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Chapter 1 


INTRODUCTION 


The main objective of the work reported in this thesis has been 
the application of fundamental concepts of continuum mechanics and 
numerical analysis for the quantitative description of some commonly 
used methods of Rapid Solidification Technology (RST) . Although 
time limitations have restricted us to the description of only one 
process in detail (namely , the Planar Flow Melt Spinning (PFMS) 
system), we believe that the information presented should allow 
the use of the same methods for the description of any other 
Rapid Solidification Processing (RSP) system. 

In this first chapter we present some background information 
about RST . After defining RSP , we comment on the merits of the 
mathematical approach to the study of RST. Then, the significant 
changes in structure and properties produced by rapid solidification 
are briefly discussed. A comment is also included regarding the 
essential process requirements for the achievement of large cooling 
and freezing rates. Furthermore, for the sake of motivation and 
completeness we briefly review Important facts concerning the actual 
and possible applications of rapidly solidified materials. We 
conclude the chapter with a summary description of what the reader 
can expect to find in succeeding chapters. 
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1.1.- Definition of RSP and the Merits of Mathematical Modeling. 

RSP is the name given to a wide array of materials processing 
operations in which the intended purpose is the production of solid 
materials directly from their melts by imposing relatively large 
cooling and freezing rates on the molten samples. Although the 
boundary between conventional casting operations and RSP is not 
clear-cut there seems to be agreement in calling RSP solidification 
processes in which the cooling rates are greater than say 10^ °C/s 
and the freezing rates larger than about 1 cm/s . 

The explicit purpose of RSP is to exploit the structure- 
properties correlation so as to be able to obtain products with 
desirable metallurgical characteristics in a reproducible fashion. 
Because of the peculiar features of RSP systems, the performance 
of suitable measurements, relatively straightforward in other 
casting systems, becomes very difficult at best. Researchers have 
resorted again and again to more unorthodox approaches in an 
attempt to reveal the fundamentals of RS processes. Mathematical 
modeling has been one of such approaches ever since the inception 
of RST. 

The mathematical modeling approach is in reality an attempt to 
make the study of RSP systems an interdisciplinary enterprise. 
When modeling, one combines the results of actual experiments with 
basic principles of continuum mechanics and concepts of numerical 
analysis to produce quantitative representations of the processes 
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under study. The explicit objectives of the approach are: 

(i) To gain an improved insight into the complex behavior of 
RSP systems, 

(ii) to establish a quantitative framework for the rational 
design and control of RST devices, and 

(iii) to demonstrate the applicability of some very basic 
principles of physics to the description of complicated, 
real-life RSP systems. 


1.2.^ The Effects of RSP on the Structure of Materials. 

The structures of rapidly solidified materials have been found 
to be markedly different to those of samples of the same materials 
processed in more conventional ways. There is ample evidence 
supporting the structure and property modifications resulting from 
RSP . Reference could be made to the review by Jones (1984) and to 
the proceedings of the various international conferences on the 
subject (see the References at the end of the report) . The interest 
generated by these discoveries is such that an international journal 
has just appeared specifically devoted to RS research. 

According to Grant (1983b), the following effects of RS are 
the main reason for the increased attention metallurgists are 
paying to RST : 

(i) Much reduced extent of segregation. The large cooling and 

freezing rates during RS do not allow time for the separation 
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and/or growth of segregated phases. As a consequence of this, 
multiphase rapidly solidified materials usually contain the 
second phases in the form of a homogeneous dispersion of very 
fine particles. 

(ii) Decreased size of microstructural features. Again , due 
to the high cooling rates involved, grains, cells and/or 
dendrites are usually much smaller than those found in samples 
produced by more conventional methods. 

(iii) The possibility of producing new phases or entirely new 
materials. It has been found that RS may prevent the formation 
of some phases commonly observed in conventionally processed 
stock. Instead, new, previously unknown metastable phases can 
appear. Moreover, the effect of RS can be so drastic as to 
entirely prevent the formation of any crystalline phases, thus 
resulting in metallic glasses. 


1.3.- Requirements for Rapid Solidification. 

There is one most important requirement that must be satisfied 
if one wishes to induce large cooling and freezing rates in molten 
samples. The requirement is the rapid formation of a thin layer of 
melt in good thermal contact with a heat absorbing medium. Rapidly 
solidified structures can also be obtained in bulk samples, however, 
if one is careful enough to eliminate all nucleating agents. Since, 
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in practice, such agents are almost always present, the requirement 
of rapidly forming thin liquid layers is really essential. Jones 
(1982) has suggested that by imposing any of the following on molten 
samples, rapid solidification can be obtained: 

(i) A high undercooling prior to solidification. Unfortunately, 
this is only possible by the avoidance of nucleating agents. 

(ii) A high withdrawal velocity of the sample through a steep 
temperature gradient. This is what is usually done during steady 
state continuous casting of rapidly solidified materials. 

(lii) A high cooling rate during solidification. This is usually 
the case during the solidification of the smaller droplets 
produced by atomization. 

From the above we can conclude that, in most cases of practical 
interest, a small sample size is required along at least one spatial 
dimension to be able to achieve the benefits of rapid solidification. 


1.4.- Properties and Applications of Rapidly Solidified Materials. 

Unexpected properties or combinations of them have been found in 
rapidly solidified materials. For example, metallic glasses tend to 
combine good ductility with high mechanical strength . Also, glasses 
with high corrosion resistance or with good catalizing properties 
have been found. Rapidly solidified stainless steels have been shown 
to be very resistant to high temperature oxidation. Useful electric 
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and magnetic properties are currently under study. The effects of RS 
on specific alloy systems are also being investigated. For example 
new alloys have been prepared by combining aluminum with unusual 
alloying elements (such as Li), resulting in improvements in thermal 
stability and mechanical properties. Property changes have also been 
reported in the cases of iron-base alloys and superalloys. These 
examples, together with the promise of more to come; have stimulated 
the current interest in RST . 

In the last few pages we have presented a summary description 
of what we believe is important background information about RST. 

At best, the information presented should enable the reader to 
understand those aspects of RST to which frequent reference is 
made in subsequent chapters without having to resort to the 
bibliography. Thus , we have presented a definition of the field 
and of the role played by mathematical models, followed by comments 
on the effects of RST on materials and the requirements for the 
achievement of high cooling and freezing rates. 

In the following chapters we take the discussion into the basic 
aspects of the modeling of rapid solidification phenomena (Ch. 2), 
then into the detailed description of the use of mathematical 
modeling techniques for the study of RSP systems (Ch. 3) . We 
conclude the report with a summary description of our results , a 
series of suggestions for further work ( Ch. 4) and the FORTRAN 
listings of the computer programs which have produced the bulk of 
the results reported ( Ch. 5) . Several appendices have also been 
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included both in order to avoid interrupting the continuity of the 
main text and to help in making the monograph more self-contained. 

We close this chapter now with a comment. Our work constitutes 
one more contribution to the rather long history of RS research 
at MIT . For the past two decades, young researchers have developed 
RS processes and have probed the characteristics of the resulting 
products. Mention should be made .of , among others , the thesis 
reports by Ruhl(1967), Strachan(1967) , Lebo(1971), Jansen(1971) , 
Domalavage (1980) , Lynch(1982), Libera (1983) , Segal(1983), Ashdown 
(1984) and Speck(1985). To them , I am sincerely indebted. Our work, 
however, has a slant in a different, relatively new direction. We 
have not performed any experiments; however, we hope to have 
demonstrated by the end of this report, that the mathematical 
modeling approach is indeed a legitimate, alternative way of looking 
at RSP problems. We believe that the development and optimization 
of RST will require a strongly interdisciplinary approach and we 
hope that our work will be regarded as an example of the way 
mathematical models can contribute to the understanding of RSP. 
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Chapter 2 


THE MODELING OF RAPID SOLIDIFICATION PHENOMENA 


In this chapter we undertake a more detailed description of 
some fundamental features of rapid solidification. Starting with 
a discussion of the effects of RS on the structure of materials 
we proceed to present summary reviews on the maximum undercoolings 
achievable in metallic melts, on the likelihood of forming metallic 
glasses and on the solidification of undercooled melts. A comment 
on the problem of morphological stability is also included. 

The attention we devote in this chapter to the undercooling 
phenomenon is due to its apparent importance in many RS systems. 
Even though our model of the PFMS process, presented in Ch. 3 , 
does not take into account undercooling effects, because the 
available evidence seems to indicate this is indeed a good 
assumption , these effects can be very significant in other RS 
systems. The prospective modeler should be aware of that, and 
should also be able to take such effects into account in his/her 
calculations, if the need arises. The information presented below 
should provide enough background to be able to follow unaided the 
current literature on the subject . 
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2.1.- More on the Effects of RSP . 

The effects produced by rapid solidification can be classified 
into two main groups: constitutional and microstructural . These 

in turn can be subdivided as follows; 

(a) Constitutional effects 

(i) Metallic glass formation. Metallic glasses can be formed 
when the rate of solidification is faster than the rate required 
for the formation of crystalline material at the solidification 
interface. Metallic glasses can be formed between: (1) metals and 
metalloids, (2) transition metal-transition metal, and (3)group II 
metals - B subgroup solutes . The kinetic conditions for metallic 
glass formation are discussed in Sec(2.3) below. 

(ii) Formation of non-equilibrium crystalline phases. It may 
well happen that the rate of solidification is faster than that 
required for the formation of the equilibrium crystalline phase. 

In this case non-equilibrium crystalline phases can appear. These 
phases can form between: (1) noble and B- subgroup metals, (2) 
B-subgroup - B-subgroup metals, and (3) transition-transition 
metals. The conditions required for the formation of non-equi- 
librium crystalline phases are described in Sec (2. 4) below. 

(iii) Solid solubility extensions. Many systems have been found 
in which RS produces significant solid solubility extension. This 
phenomenon was indeed the first indication of the effects of RS. 
In Sec (2. 5) we discuss solute redistribution during RSP, 
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(b) Microstructural effects 

(i) Grain structure and size. Equiaxed grains, cells and 
dendrites have all been observed in rapidly solidified specimens. 
However, in all cases the sizes involved are much smaller than 
those that can be obtained by conventional processes. In Sec (2. 2) 
we review the thermal conditions required for the production of 
microcrystalline structures from the melt and in Sec (2. 6) we 
briefly describe the main features of the problem of morphological 
stability during RS . 

(ii) Formation of lattice defects. Few heavily dislocated 
specimens have been found as a result of RS. Few twins and some 
stacking faults have also been observed. On the other hand, 
considerable vacancy supersaturations have been found. 

The presence of some or all of these effects in rapidly solidified 
materials accounts for their unexpected properties. In the sequel 
we present a summary description of the phenomena which produce 
such effects. 


2.2.~ On the Maximum Attainable Undercooling During RSP. 


Since the solidification structures of heavily undercooled 
samples are usually composed of very finely grained material it 
has been natural to ask about the maximum possible undercooling 
a given sample can achieve prior to solidification. The problem 
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has been studied by Hirth(1978) using concepts from nucleation 
theory. It is well known that, given sufficient undercooling prior 
to solidification, nucleation and growth can proceed to complete 
solidification without the solid-liquid interface ever reaching 
the solidus temperature . A sample in this condition is said to be 
hypercooled. The condition for hypercooling can be derived from 
thermodynamics and it is 

=p « T L - V + (T s - t l» - L (1) 

where all the symbols are defined in the List of Symbols. We note 
that Eqn(l) is only valid for the situation in which the temperature 
is uniform across the sample (Newtonian cooling), however , despite 
this limitation it constitutes an useful estimate. It should also be 
mentioned that under hypercooled conditions, the rate controlling 
step of the solidification process are the atomic jumps across the 
solid-liquid interface . 

The maximum attainable undercooling is, obviously, the one 
corresponding to homogeneous nucleation. An index of merit can then 
be defined for the achievement of nucleation control. If the forma- 
tion of one nucleus per sample constitutes a suitable definition of 
the critical condition for the production of homogeneous micro- 
crystalline structures, the use of homogeneous nucleation theory 
leads to the following expression for the maximum attainable under- 
cooling, 
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i6 tt <r 3 o 2 


T - T 
L N 


1/2 

) 


3 kg T n L Z ln( 10 3 d 3 (r*/a) Z (a/0)D 1 (T L - T N )/T ) 


( 2 ) 

Since there are only small differences in the nucleation 
behavior of single and multicomponent systems, Eqn(2) can be used to 
estimate maximum undercoolings in both cases. For a given material 
and sample size , Eqn(l) can be used to estimate the minimum under- 
cooling required for hypercooling. Afterwards, Eqn(2) can be used 
to compute the minimum cooling rate required for the production 
of homogeneous microcrystalline structures. 

Hirth has been able to obtain reasonable estimates of critical 
cooling rates for various alloy systems. His final advice is , 
however, that more attention be paid to the role of impurities and 
other heterogeneities in promoting nucleation, thus preventing the 
attainment of the maximum undercooling. 


2.3.- Metallic Glasses. 

When the cooling rates during RSP are sufficiently large , 
metallic glasses can be formed instead of the usual crystalline 
phases. The formulation of the critical conditions for metallic 
glass formation has been worked out by Uhlmann(1983) . 
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From the Johnson-Mehl-Avrami theory of phase transformation 
kinetics, the relationship between the fraction of new phase formed. 


X and the time t is 
c 


X = (1/3) 7f I R 3 t 4 

C V 


( 1 ) 


The nucleation rate 1^ can be estimated from homogeneous nucleation 
theory and is 

1.024 

I v = N° 9 exp ( ) (2) 

(T/T f ) 3 ((T f -T) /T f ) 2 


The growth rate can be estimated from crystal growth kinetics to be 


R = 0.2 ( (T f - T)/T f ) V> a( 1 - exp(- 


L((T, - T) fl c ) 


) ) 


k N T 
B A 


(3) 

Finally , the jump frequency V , can be related to the viscosity 
of the melt through the Stokes-Einstein relationship, i.e. 



C4) 
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Equations (1)^(4) above can now be used to construct temperature- 
time-transformation diagrams for the prediction of the critical 
conditions for metallic glass formation. The procedure used for the 
construction of such diagrams is as follows: 

(i) First , an arbitrarily small fraction crystallized (say 

X* = 10 ^ ) is selected and a temperature (below the solidus) 
c 

is chosen. 

(ii) From Eqns(l)-(4) the time required for fraction X* to 

form is calculated . The pair T , t is plotted on a scale of 

T vs. log t . 

(iii) Steps (i) and (ii) above are repeated for different T f s 

(but the same X* ) to obtain the complete T-T-T diagram, 
c 

Figure (1) shows the result of one such calculation for various 
metallic systems. Based on the same ideas, Uhlmann estimated the 
maximum thickness of sample that can be transformed into glass by 
rapid cooling to be 


1/2 

y g = ( * t c ) (5) 

where t is shown in Fig.(l) for the case of AuGeSi. 
c 

Other factors that have been found to influence the formation 
of metallic glasses are , heterogeneous nucleants and nucleation 
transients . For a discussion of these factors the reader shoul see 
the paper by Uhlmann mentioned above. 
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T f (Ni) 



log t 


Fig(2.3.1) Temperature-time-transformation (T-T-T) diagrams 
for the crystallization of several metals from their undercooled 
melts. Here X* = 10 ^ . From Davies (1976) . 
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2 .4.’- Solidification of Undercooled Melts. 

When a melt sample which has been undercooled starts to solidify, 
the latent heat of solidification is released very rapidly at the 
solid-liquid interface. The sudden release of this latent heat is 
so fast that the outer surface of the sample may well be unable to 
dissipate this energy. The liberated heat has to be retained in the 
melt thus producing the phenomenon known as recalescence . During 
recalescence , the temperature of the solidification interface 
rises quickly and can even reach the equilibrium value. 

The solution of the coupled heat transfer - crystal growth 
problem is not a simple task. Solution procedures usually start 
by assuming a particular expression for the crystal growth rate 
as a function of the interface temperature. The heat transfer 
problem is handled essentially in the same way as the usual Stefan 
problem (see Sec(3A.2)). However , instead of the fixed , equilibrium 
freezing temperature found in the classical Stefan problem, here 
the interface temperature is variable. 

Two growth rate - interface temperature relationships have been 
the most popular, namely , the exponential law 


R = R ( 1 - exp( L(T - T,)/T T, ) ) 
o t r 


(i) 


and the linear law 


R = R^ (T - T f ) 


( 2 ) 
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The solidification problem is described by the equations of the 
classical Stefan problem in the weak (enthalpy) formulation (see 
Sec(3A.2)) , 


t = div( K grad T ) + r^ (3a) 

E = f (T) (3b) 

However, instead of assuming the solid-liquid interface temperature 
to be given by equilibrium considerations, it is assumed to be 
dependent on the crystal growth rate according to either Eqn(l) or 
Eqn(2) above . 

The mathematical problem represented by either Eqns(l) and (3) 
or (2) and (3) must be solved to obtain the temperature field inside 
the sample, the freezing rate and the interface temperature. Various 
methods have been proposed for the solution of this problem. 

Boswell (1979) used a front-tracking technique (Sec (3A. 2) ) to 
solve Eqns(l) and (3) for the case of a pure metal solidifying 
against a metal chill. An iterative method was used to compute the 
interface temperature. Although no details were given, he presented 
a plot describing the effects of the heat transfer coefficient, of 
the splat-chill interface temperature at the start of freezing and 
of the materials properties. 
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Levi and Mehrabian(1982) performed a detailed analysis of the 
rapid solidification phenomena taking place during the freezing of 
undercooled pure metal droplets. They solved Eqns(l) and (3) and 
(2) and (3) in a suitable coordinate system using appropriate 
boundary conditions. The most important parameters in their cal- 
culations were the droplet diameter, the heat transfer coefficient, 
the droplet surface temperature at the onset of nucleation, and the 
kinetic growth coefficients and in Eqns(l) and (2). 

Two heat transfer models were used by Levi and Mehrabian. The 
simplest one assumed negligible temperature gradients inside the 
droplet (Newtonian model). In the other model , this restriction 
was relaxed. An implicit finite difference method was used to 
solve the equations for this latter case. Although somewhat similar 
results were obtained from both models, the one based on non- 
Newtonian cooling provided more detailed information. The results 
of their calculations were conveniently summarized in the form of 
dimensionless enthalpy-temperature curves. The cooling and freezing 
paths of individual droplets can be easily followed in such diagrams. 
One example of such plots is included in Fig(l). We will describe 
now how to read this diagram. 

At the lower limit of slow cooling rates, the freezing path 
followed by the droplets is close to the equilibrium freezing path. 

In this case, when the droplet reaches the equilibrium freezing 
temperature T^ , solidification at constant temperature starts 
and continues until the entire latent heat has been withdrawn ( 
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T (-) 

Fig(2.4.1).- Dimensionless enthalpy-temperature diagram showing 
the various possible freezing paths of undercooled melts. From 
Mehrabian(1982) . 
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and the dimensionless enthalpy is equal to zero). This solidification 
mode is called isothermal for obvious reasons. 

At high cooling rates the other extreme possibility appears. 

Namely , during cooling, an amount of energy at least equal to the 
latent heat of fusion is withdrawn without nucleation taking place. 
Once this is done the system can solidify without having to extract 
any more heat from the sample. This is the so called isenthalpic 
(or adiabatic) solidification mode. 

A much more common occurrence is the intermediate case where 
solidification starts below the equilibrium freezing point but 
before the entire latent heat of solidification has been released. 
Since nucleation is accompanied by the liberation of a certain 
amount of heat, the droplet temperature will tend to rise until 
the equilibrium melting point is almost reached. Solidification 
can then proceed according to the isothermal mode. This self- 
heating process is known as recalescence. 

Two distinct solidification regimes can thus be observed in 
this the more general case. First, during recalescence, the 
solidification interface moves rapidly into the liquid. The latent 
heat is released so quickly that the external cooling is unable 
to extract it thus resulting in the heating up of the droplet. 
However, afterwards, when the droplet temperature has reached a 
value close to the equilibrium melting point, the subsequent 
freezing depends mainly on the rate of heat extraction through the 
outer droplet surface . 
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As expected, the microstructures of the products formed during the 
two freezing stages in the general case, are markedly different. A 
detailed review of the subject , including many photomicrographs, has 
prepared by Mehrabian. 

More recently, Crowley(1984) studied the process of pulsed laser 
annealing. The large heating and cooling rates produced by this 
process induce undercooling. She proposed a fixed domain method with 
partial front tracking for the solution of this problem. The 
formulation is again identical to that of the conventional Stefan 
problem except for the allowance of a variable solidification 
interface temperature. Crowley produced a consistent and stable 
algorithm free of oscillations. Her technique certainly warrants 
attention from people studying the solidification of undercooled 
melts . 

In a related study, Dantzig and Davis(1978) used the same basic 
set of equations in combination with alternative mathematical 
methods (matched asymptotic expansions with embedding) to analyze 
the conditions for non-equilibrium phase formation during RSP. 

They introduced the concept of the delay time as the time that 
must elapse between the attainment of the equilibrium melting 
point and the moment when the melt transforms into the equilibrium 
product. By comparing the delay time with the times required for 
non-equilibrium products to form, they derived a criterion for the 
formation of the latter. The basis for comparison was the difference 
in rates of the process of interfacial attachment and the process of 
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heat conduction. The exponential law for crystal growth (Eqn(l)), 
was solved simultaneously with the equations of the classical 
Stefan problem modified by the presence of the delay time in the 
Stefan condition. 

The conclusion reached by Dantzig and Davis was that, if the 
kinetic processes of atomic attachment at the solidification 
interface are slow compared to the cooling rate, the expected 
equilibrium crystalline phase may never form. Instead, a super- 
cooled layer of fluid will grow from the chill until it reaches 
macroscopic dimensions. The critical delay time separating 
equilibrium from non-equilibrium products was calculated to be 

= - < V ^ ) (*> 

Equation (4) shows the expected result, that low nucleation 
temperatures and high cooling rates facilitate the formation of 
non-equilibrium products. 

Clyne(1984) has presented a review of the numerical treatment 
of RS processes in which undercooling is an important consideration 
and his paper can be consulted for additional information. 


2.5.- Solute Redistribution During RS . 


The phenomenon of solute redistribution during RS is still the 
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subject of active research. Although many aspects of it are still 
obscure, important insight was gained from the model proposed by 
Kattamis (see e.g. Flemings (1981) ) . This model suggests that the 
freezing of undercooled alloys takes place according to the 
following three stages: 

(i) Recalescence up to the solidus temperature , 

(ii) recalescence from T up to the maximum recalescence 

temperature T* , and 

(iii) cooling from T* . 

In the model it is also assumed that the diffusion of solute 
is negligible during stage (i) but not during (ii). The segregation 
during stage (iii) is described by the Scheil equation. Next we 
present a brief summary of the equations of this model. 

For stage (i) above, a heat balance can be written as 


d f /d T = C /L 
s p 


( 1 ) 


so that the fraction solidified once the solidus temperature is 
reached during recalescence, f* , is 


f s ’ (C p /L) (T S - V 


( 2 ) . 


Since diffusion is neglected during this stage, the solute 

concentration in the solid forming between T and T is 

NS 


C = = C* 

s • s 


( 3 ) 
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For stage (ii), a solute balance can be written as 


d f /d c. = (1 - f )/(C - C*) (4) 

S -L S X. S 

Moreover, from the phase diagram, the slope of the solidus curve 
is given as 


m = d T/d C* (5) 

o s 

The combination of Eqns(4) and (5) leads to 

d f /d C* = (C /L) m (6) 

s s p S 

which can be integrated between f* and f to give 

s s 

C* - C„ + (t/. s c p )(f s - f*) (7) 

Substituting now Eqn(7) into Eqn(4) and assuming that f = 

i a. 

f and that (1 - f ) = constant , leads to 
s s 

± ^ 

(C x - )(f s -l) " (f s - f s ) ( L/(2 m s C )) - 0 (8) 

Now, since C* = k when T = T* , form the assumption 
of local equilibrium at the solidification interface, the combina- 
tion of Eqns(7) and (8) allows us to compute the fraction 
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solidified when the maximum recalescence temperature is reached. 



Finally , stage (iii) is assumed to take place according to the 

Scheil equation modified by the fact that the initial state is 

given by f = f** . The resulting expression is 

s s 


k - 1 

c* = c* 1 ^ i - ((f s - fJVa - fg 1 )) ) (9) 

Where C*^ is the solute concentration in the solid side of the 
s 

solidification interface corresponding to the maximum recalescence 

temperature T* . Because of the assumed negligible diffusion in 

the solid, all the interfacial concentrations mentioned above are 

basically equal to the final concentrations inside the solid, l.e. 

C = C* . 
s s 

The model described above provided the first quantitative 
explanation for the frequently observed solute rich cores of 
dendrites found in samples produced by solidification of under- 
cooled melts. The model has been refined and alternative stages 
have been proposed. The thesis by Chu(1983) contains a detailed 
description of the state of the art in this area and it should 
be consulted for further information. 
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2 . 6 .- 


Morphological Stability During RS 


Since material properties are the main concern of the 
metallurgist and because these properties are strongly related 
to the microstructure, the prediction of the relationship between 
the process parameters and the resulting microstructure has long 
been an important consideration. This has also been the case in 
RS research. The main question to be answered is if the solid- 
liquid interface will grow in a planar fashion without micro- 
segregation or will it break up into cells or dendrites, resulting 
in segregated, multiphase structures. 

The principle of constitutional supercooling provides a useful 
guide to ascertain the growth conditions resulting in solidification 
interface shape instability during alloy solidification. However, 
research on RS has shown that the constitutional supercooling 
principle produces entirely erroneous predictions in the extreme 
case of large freezing rates. This deficiency has been removed by 
the introduction of the theory of morphological stability (see e.g. 
Coriell and Sekerka(1980) or Cahn et al.(1980)). 

The theory of morphological stability is based on a kinetic 
analysis of the spatial and temporal behavior of perturbations 
formed on the solidification interface. The starting point for 
the analysis is the governing equations for heat and mass transfer 
(see Sec(3A.l)). A perturbation -type linear stability analysis is 
employed to derive the conditions for stability. The simplest 
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model adopted for study is the directional solidification of a 
binary alloy under a constant growth velocity and under local 
equilibrium conditions at the solid-liquid interface. 

The governing equations for heat and mass transfer must account 
for the latent heat released during solidification and thus they 
are basically the same as those of the Stefan problem for alloy 
solidification (Sec(3A.2)), except for the incorporation of 
interface curvature effects in the boundary conditions. First, the 
temperature field , the concentration field and the interface shape 
are written as the product of a constant term and a perturbed part, 
i.e. 


T 


C 


F 


T exp( At + 1( CO x + 

o X 

<O y y) ) 

(l) 

C exp( A t + i( CO x + 

o X 

CJ y y) ) 

(2) 

F q exp( A t + i( ^0 x x + 

<^> y y) ) 

(3) 


From the form of these expressions, it can be readily seen that the 

interface will become unstable whenever the real part of A i s 

positive for any real values of CO and CO . It is possible 

x y 

to derive an equation for the quantity A as a function of the 
process parameters and the material properties by simply substi- 
tuting Eqn(l)-(3) into the original governing equations. The use 
of appropriate boundary conditions leads finally to the desired 
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expression which will not be quoted here but can be found in the 
references given. The main feature of this equation is , however, 
that it is composed of three terms: a term involving thermal 
effects, another involving surface tension effects and the last one 
involving concentration effects. From the form of the equation 
it is seen that the thermal and surface tension effects tend to 
dampen the interface shape perturbations and are thus stabilizing. 
The concentration term , on the other hand, has always a de- 
stabilizing effect. When this last term is sufficiently large, the 
interface may become unstable and perturbations will grow. 

The stability equation mentioned above can be simplified con- 
siderably if due account is taken of the extremely large thermal 
'diffusivities of metallic melts by assuming it equal to infinity. 
Under this assumption the stability equation becomes. 


2 Gj + K L O - ( R s + V G c S(A g ,k) (4a) 

where 

Gj_ - (d T/d x) 1 (4b) 

G = R C (k - 1)/D. k (4c) 

c 1 

S(A ,k) = 1 + (A /4k) (1 - r2 + 2kr 2 ) - (3 k 1 ' 1 r/2) (4d) 

s s s 

A s = k T f ( G/ j> L) R 2 /D 2 G c (4e) 
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and the quantity r must be obtained by solving 


r 3 + (2k - 1) r - (2k/A 1/2 ) = 0 (4f) 

s 

Equations (4a) — (4 f ) provide a suitable representation of the 
stability behavior of metallic melts for a wide variety of process 
conditions. Two limiting cases exist, namely, for small growth 
velocities , the constitutional supercooling criterion is adequate 
and the stability limit is thus given by (Flemings(1974) ) , 


G /R = - m C* (1 - k)/k D 

1 L S 1 


(5) 


On the other hand, for large growth velocities, the absolute 
stability limit (obtained by making A g = 1 in Eqn(4)), is a 

good approximation, i.e. 

G c /R = k T f (<T/ j> L) m L R/D 1 (6) 

Equation (6) indicates that much greater stability can be 
expected at high freezing rates than that predicted from the 
constitutional supercooling criterion. At least three factors 
can explain this behavior. First, the capillary forces have a 
strongly stabilizing effect, particularly on the short wavelength 
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perturbations characteristic of high growth rates. Furthermore, 
the deviations from local equilibrium at the solidification inter- 
face and the peculiar temperature gradients resulting from the 
freezing of undercooled melts also have a stabilizing effect. 

A convenient way of presenting the results of the calculation* 
performed using Eqn(4) or (5) and (6) is by plotting the bulk 
liquid solute concentration against the growth velocity. The pairs 
of values of these quantities corresponding to the critical 
condition display the limit of stability. One such plot, for the 
case of the Al-Cu system is presented in Fig(l). 

The theory of morphological stability has been extended to 
deal with other effects such as interfacial anisotropy, melt 
undercooling, nonlinearities , and deviations from local 
equilibrium at the solidification interface. For this latter case, 
a corresponding stability equation has been derived using ideas 
very similar to those described above. The analysis has provided 
useful insight about the important phenomena of solute partitio- 
ning and trapping during RS . The references should be consulted 
for details. 

In this chapter we have reviewed several topics concerning the 
mathematical representation of rapid solidification phenomena. 

As can be inferred from the discussions on undercooling, metallic 
glasses and the freezing of undercooled melts, the fundamental 
processes of nucleation and growth played an important part in the 
description of such systems. However, when kinetic considerations 
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Fig (2 f 6 . 1) , - Interface stability diagram for the directional 

A 

solidification of Al-Cu alloys. Here = 2* 10 K/m . From 

Coriell and Sekerka(1980) . 
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of this type are tried for the description of the more complex 
RSP systems found in practice, the mathematical problem becomes 
very difficult. The non-constant growth rates, the poorly defined 
boundary conditions and the existence of undetermined computational 
domains, all contribute to the difficulties. 

A very useful simplification is obtained when the kinetic 
processes taking place at the solidification interface are assumed 
to be so fast that they can be safely disregarded as the rate 
controlling step of the overall process. Under these circumstances, 
the macroscopic transport processes control the overall performance 
of the system. In the following chapter we present some solutions 
to the mathematical problem resulting from the neglect of the 
atomic kinetic processes at the solidification interface. Only 
one system ( the PFMS ) is dealt with in all detail and summary 
comments are included for a few others. 


Chapter 3 


THE MATHEMATICAL MODELING OF RAPID SOLIDIFICATION PROCESSING 


As mentioned at the end of the previous chapter, the solution 
of problems in RST can be facilitated if the atomic kinetic 
processes taking place at the solidification interface are 
assumed to be so fast that they can safely be disregarded as the 
controlling mechanism for the overall process. This is equivalent 
to assuming that the rate controlling processes are of a macro- 
scopic nature. It is indeed fortunate that the assumption of 
infinitely fast interfacial processes is justified for substances 
constituted by small molecules (such as metals) in many cases of 
practical interest . 

In this chapter we undertake the task of simulating mathemati- 
cally the behavior of some important RSP systems using the 
assumption of infinitely fast interfacial processes. For the sake 
of organization, in the first section we attempt a classification 
of RSP systems which is capable of including every existing (and 
non-existing) rapid solidification technique. We then proceed to 
the formulation of the simplest macroscopic heat transfer models, 
based on the assumption of Newtonian cooling conditions. These 
simple models can be very useful to obtain first order estimates 
of cooling and freezing rates in actual RSP configurations. 
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Although the simple heat transfer models have been widely used, the 
interpretation of the subtle variations found from process to process 
requires models of greater accuracy. 

In Sec (3. 3) we describe the somewhat more sophisticated models 
resulting from the elimination of the assumption of Newtonian 
cooling. Since the details of these models are highly system- 
specific, only one RSP system is dealt with in all detail while 
the basic ideas required for the formulation of the models of other 
important systems are the subject of much briefer presentations. 

We focus our attention on the Planar Flow Melt Spinning System 
(PFMS) and present enough detail, that the extension of our methods 
to other RSP systems should be relatively straightforward . 

To facilitate the reading we have decided to separate background 
information from that pertaining specifically to the modeling of 
RSP systems. However, for completeness , the background informa- 
tion has been put in appendices at the end . 


3.1.- Rapid Solidification Processing Systems. 

A large variety of devices have been constructed and used for 
the production of rapidly solidified materials. Most of them, 
however, have been designed having in mind the main requirement 
for obtaining high cooling and freezing rates, namely, the 
existence of a small section in at least one spatial direction. 
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Jones(1982) has proposed a classification of RSP systems based 
on some key features of the various processes. He considers RSP 
systems to be divided into: (i) Spray methods (involving the 

complete disruption of the continuity of the melt), (ii) chill 
methods (where the melt is thinned instead of being disrupted), 
and (iii) weld methods (where a high energy beam melts the surface 
of a bulk object). We believe that the classification presented 
in Table (1) below, which is based on Jones' , is more comprehensive 
and it is the one we will use in our discussion. 


Actual, representative examples of all the categories listed 
in Table (1) can be found in Table (2) together with references 
where the actual devices are described. To aid in the reading of 
Table (2) , Fig(l) shows schematically some of the most important 
processes included in this table. Interestingly enough , most of 
the seemingly entirely different processes included in Table (2) 
have important features in common. The basic physical phenomena 
involved with the performance of RSP systems are described in 
Table (3) . A glance at this table readily shows that the most 
important physical processes taking place during RSP operations 
are: (i) The fluid flow phenomena associated with the spreading, 

squeezing, thinning and breaking up of molten metal samples, and 
(ii) the energy transfer processes governing the cooling and the 
solidification of such samples. 

It should be noted that even though very much the same physical 
processes are at work in all RSP systems, subtle differences in 


design can produce significant changes in the characteristics of 
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Table (3,1.1).- Rapid Solidification Processing Systems 


I) Melt Fragmentation Processes 

a) Fragmentation produced by moving solids 

b) Fragmentation produced by moving fluids 

II) Splatting Processes 

a) Splatting to produce particulate material 

b) Splatting to produce continuous material 

III) Direct Quenching Without Fragmenting or Splatting 

IV) Melting and Quenching of Thin Surface Layers 

V) Liquid Dynamic Compaction and Spray Deposition 
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Table (3.1.2).- Some Examples of RSP Systems 


la) Melt Fragmentation Produced by Moving Solids 


1) Rotating Cup or Dish 

2) Rotating Perforated Cup 

3) Rotating Electrode Process 

4) Impact Disintegration 

5) Vibrating Electrode 

6) Melt Drop Technique 

7) Twin Roll Technique 

8) Single Roll w/Serrated Surface 

9) Melt Extraction w/Serrated Wheel 

10) Single Roll Atomization 


Glickstein et al(1978) 
Daugherty (1964) 

Champagne & Angers (1984) 
Schmitt(1979) 

Ruthardt & Lierke(1981) 
Aldinger et al(1977) 
Singer et al(1980) 
Carbonara et al(1982) 
Pond et al(1976) 
Narasimha & Sekhar(1984) 


lb) Melt Fragmentation Produced by Moving Fluids 


1) Water Atomization 

2) Subsonic Gas Atomization 

3) Ultrasonic Gas Atomization 


Tallmadge(1978) 
Beddow(1978) 
Grant (1983a) 
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Table (3.1.2).- (contd.) 


Ila) Splatting for Particulates 

1) Gun-Ski jump Device 

2) Piston and Anvil Device 

3) Injection Chill Mold 

4) Isolated Droplets on Chill 

5) Rotating Impactor 


Duwez & Willens(1963) 
Strachan(1967) 

Hinesley & Morris (1970) 
Madejski(1976) 

Predel (1978) 


lib) Splatting for Ribbon or Sheet 

1) Chill Block Melt Spinning 

2) Centrifugal Melt Spinning 

3) Planar Flow Melt Spinning 

4) Melt Drag 

5) Twin Roll Quenching 

6) Melt Extraction 


Liebermann & Graham (197 6) 
Chen & Miller(1976) 
Fiedler et al(1984) 

Hubert et al(1973) 

Murty & Adler(1982) 
Robertson et al(1978) 
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Table (3.1.2).- (contd.) 


III) Direct Quenching Without Fragmenting or Splatting 

1) Melt Extrusion Shepelsky & Zhilkin(1968) 

2) Taylor Wire Process Manfre et al(1974) 

3) Free Flight Melt Spinning Kavesh(1976) 

IV) Melting and Quenching of Thin Surface Layers 

1) Laser Processing Breinan & Rear (1983) 

2) Electron Beam Processing Mawella(1984) 

V) Liquid Dynamic Compaction and Spray Deposition 

1) Liquid Dynamic Compaction Singer & Evans(1983) 

2) Plasma Deposition Apelian et al(1983) 
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Table (3.1.3).- Basic Physical Processes During RSP 


la) Melt Fragmentation Produced by Moving Solids 

1) Fluid Flow Phenomena 

a) Impact and Spreading of Melt on Moving Solid 

b) Melt Thinning and Acceleration 

c) Melt Fragmentation Proper 

i) Direct Drop Formation 

ii) Ligament Formation 

iii) Film Formation 

d) Bursting of Melt by Impactor 

e) Capillary Wave Atomization 

f) Cavitation Inside Melt 

g) Shearing of Melt by Serrated Disk 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 
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Table (3.1.3).- (contd.) 


Ib) Melt Fragmentation Produced by Moving Fluids 

1) Fluid Flow Phenomena 

a) Melt Thinning 

b) Growth of Disturbances on Melt Surface 

c) Formation and Tearing of Ligaments from Melt Sheet 

d) Growth of Disturbances on Surface of Ligaments 

e) Formation and Separation of Droplets 

f) Droplet Breakup 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 

Ila) Splatting for Particulates 

1) Fluid Flow Phenomena 

a) Impact and Spreading of Melt on Substrate 

b) Squeezing of Melt between Two Substrates 
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Table (3.1.3).- (contd.) 


2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 

Tib) Splatting for Ribbon or Sheet 

1) Fluid Flow Phenomena 

a) Ejection of Melt from Nozzle 

b) Impact, Adhesion and Spreading of Melt on Moving Chill 

c) Impact, Adhesion , Spreading and Squeezing of Melt 
between Nozzle and Moving Chill 

d) Squeezing of Melt betwee Two Moving Chills 

e) Capillary Flows 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 
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Table (3.1.3).- (contd.) 


III) Direct Quenching Without Fragmenting or Splatting 

1) Fluid Flow Phenomena 

a) Stabilization of Liquid Metal Jet 

b) Velocity Relaxation in Melt Jet 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 

IV) Melting and Quenching of Thin Surface Layers 

1) Fluid Flow Phenomena 

a) Motion on Free Surfaces 

b) Surface Tension Driven Flows 

c) Natural and Forced Convection 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 
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Table (3.1.3).- (contd.) 


V) Liquid Dynamic Compaction and Spray Deposition 

1) Fluid Flow Phenomena 

a) Impingement, Spreading and Mixing of Falling Droplets on 
either Shallow Melt Pools, Mushy Surfaces or Solid Droplets 

2) Heat Transfer Phenomena 

a) Cooling 

b) Freezing 
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the products of processing. The varying degrees of interaction 
between the fluid flow and the heat transfer phenomena in the 
various processes account for the observed differences in process 
performance. For example, although a molten jet is thinned during 
both gas atomization and melt spinning, complete breakup to form 
powder is the final objective in the first case, while the forma- 
tion of a continuous ribbon is desired in the latter. It is, thus, 
the interplay between spreading and thinning rates on the one hand 
and cooling and freezing rates on the other that accounts for the 
wide variety of existing RSP systems . 

As expected, the different techniques will produce rapidly 
solidified products with structures (and properties) peculiar to 
them and thus, widely different microstructures may be found in 
samples of the same material produced by different techniques. 

This complexity makes necessary a case by case study of the 
various processes. Fortunately, despite the idiosyncracies of the 
different RS techniques, the same fundamental principles of 
continuum mechanics are applicable to all of them. This fact 
provides the unifying feature for the mathematical modeling of 
RSP systems. 


3.2.- Mathematical Models for RSP Systems. Newtonian Cooling. 


Mathematical models based on heat transfer considerations have 
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long been used to estimate the cooling and freezing rates obtained 
during RSP . Because of their inherent simplicity, lumped parameter 
models were almost always invariably adopted. The main assumption 
involved in all of these early models was the neglect of temperature 
differences across the sample thickness , i.e. Newtonian cooling 
conditions. Mathematical models of heat transfer and solidifi- 
cation based on the assumption of Newtonian cooling always lead 
to ordinary differential equations which are relatively easy to 
solve. 

In this section we describe the formulation and the solution 
of mathematical models of RSP systems based on the assumption 
of Newtonian cooling. Because of their simplicity, the models can 
be very general. Furthermore, to avoid the drudgery of hand 
calculating the cooling and solidification rates we have included 
(in Ch. 5) a computer program capable of doing all the necessary 
computations . 

In the description which follows we first present the for- 
mulation for the processes resulting in separated particles and 
then go on to describe the formulation for those processes 
resulting in ribbon, sheet or wire. 

a) Lumped parameter models for discrete splats 

Let A and V be , respectively, the splat surface in contact 
with the heat sink and the total volume of the sample (Fig(l)). An 
overall heat balance for the splat is composed of three separate 
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Fig(3 .2.1) Schematic representation of RSP systems used for 
heat transfer calculations according to the lumped parameter model 
(a) Sphere , (b) cylinder, and (c) slab. 


8 




stages: 

(i) Cooling of the melt down to the melting point 


J> V C p dT/dt = - h (T - ) A 


( 1 ), 


(ii) solidification of the sample at constant temperature T 


f 


vy> L df g /dt = h (T f - T* ) A (2), 

and 

(lii) cooling of the solidified sample 


J> V C p dT/dt = - h (T - ) A (3). 

Integrating Eqns(l)-(3) between suitable limits produces the 
following expressions for the temperature and the cooling and 
freezing rate; 

For stage (i) 

T = (T p ~ T oo > ex P< ~ h A t/V J> C p ) + (4a) 

dT/dt = - (T p - T — ) (h A/V p C p ) exp ( - h A t/V C ) (4b) 

For stage (ii) 

f s - h A (T f - T„ )<t - t ss )!VJ> L (5a) 
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h A (T f - T* )/V £ L 


(5b) 


df / dt 
s 


And for stage (iii) 

T - (T f - ) exp( - h A (t - t eg ) /V p C p > + T * < 6a > 


dT/dt = - (T f - T )(h A/V p C p ) exp(- h A (t - t es >/V/> C p ) 

(6b) 

In these expressions t is the time for the start of freezing 

s s 

and t is the time for the end of solidification, 
es 

If we write Z for the radius of the sphere or of the 
cylinder or for the half -thickness of the slab in Fig(l), the 
following relationships hold. 


A/V = 

3/Z 

for 

the sphere 

(7a) 

A/V = 

2/Z 

for 

the cylinder 

(7b) 

A/V = 

1/Z 

for 

the slab 

(7c) 


Moreover , the solidified thickness at any given time can be 
obtained from the fraction solidified as follows, 

m 

z = Z (1 - f g ) (8) 
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where m takes the values of 1/3 , 1/2 , and 1 , respective- 
ly, for the sphere, the cylinder and the slab. We note that all 
these relationships were derived for the situation in which the 
heat extraction takes place through all sides of the sample. 

However , the heat lost through the ends of the cylinder or 
through the edges of the slab has been neglected. 

Equations (4)-(8) can be used to estimate cooling and freezing 
rates in a wide variety of RSP systems. The FORTRAN program 
RSPNN presented in Sec (5.1) below has been designed to perform 
these calculations. 

b) Lumped parameter models for continuous processes 

Let H(x) be the local melt thickness (see Fig(2)). In this 
case we perform the overall heat balance on volume elements of 
size Ax along the downstream direction. These volume elements 
are assumed to be moving in a plug flow fashion. Proceeding as 
before, after integration , the following expressions for the 
temperature and the fraction solidified can be obtained, 

T x + 4x - < T x - T - > exp( - h A Ax/V > > C p ? x > + T - 

(9) 

and 

f s - f + h A (T f - T.. ) &x/Vf L V x (10) 

x + Ax x 
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(a) Sheet cooled from one side, and (b) from two sides 


■ 


where the subindex x denotes the location along the downstream 

direction corresponding to a given T or a given f . 

s 

From geometrical considerations again , the quantity A/V is 
given by 


A/V = 4/H for the cylinder 

A/V = 2/H for the sheet cooled from two sides 

and 


A/V - 1/H for the sheet cooled from one side. 

Moreover, the relationship between the actual solidified thickness 
and the local fraction solidified is , 


1/2 

y_ = (H/2)(l - (1 - f ) ) for the cylinder 

£> S 


y = (H/2) f for the sheet cooled from two sides 

s s 


and 


y ~ H f for the sheet cooled from one side, 
s s 


Equations (9) and (10) must be applied repeatedly, marching 
forward along the downstream direction to obtain cooling and freezing 
profiles. 

The single most important adjustable parameter in the above 
equations (as well as in those to be presented in Sec(3.3) below), 
is the heat transfer coefficient. The value of this coefficient 
directly determines the cooling and freezing rates of the rapidly 
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solidified sample. Since these rates are directly related to the 
structure of the material, the availability of realistic values of 
h is a question of great importance. The use of heat transfer 
coefficients to describe the complex heat transfer phenomena which 
take place at interfaces between phases is justified as long as 
there are no more rigurous means of describing such processes. 

For the convenience of the reader who wishes to perform heat 
transfer calculations similar to the ones reported in this thesis, 
and also for comparative purposes, we have compiled in Table (1) 
a list of suggested values of the heat transfer coefficients for 
a wide variety of casting/solidification processes. The sources 
have also been included. 

In the past, thermal and microstructural measurements have 
been employed for the estimation of h .We now suggest the use 
of the coefficients in Table (1) to calculate thermal responses 
and the resulting microstructures . 


3,3.- Mathematical Models for RSP Systems. Non-Newtonian Cooling. 

In some cases , the assumption of Newtonian cooling conditions 
can be grossly inadequate. It may be that the temperature gradients 
across the splat just cannot be neglected. This is particularly 
true for those RSP systems whose performance strongly depends 
on the existence of large temperature gradients across the splat. 
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Table (3.2.1).- Heat Transfer Coefficients Typical of Casting- 


Solidification Processing Operations. 


System 

Sample Size 

h 

Source 


( y*m) 

(cal/cm^s°C) 


Rotating Dish 



Glickstein 

Atomization of 

20 - 500 

0.2 - 700 

et al (1978) 

IN-100 (P&W) 




Gun-Ski jump 
Splat of A1 on Ni 

0.1 - 5 

2.7 - 6.8 

Predecki 
et al (1965) 

Piston & Anvil 
Splat of A1 on Fe 

76 

0.4 - 5 

Harbour 
et al (1969) 

Die Casting 
A1 on Steel 

1600 

1.9 

Mehrabian 

(1982) 

Metal Splat on 
Metal Substrate 

not given 

2.4 - 24 

»» II 

Atomization of A1 




(Radiation Cooling) 

100 

0.0013 

Jones(1982) 

Atomization of A1 




(Convection Cooling) 

100 

0.24 - 2.4 

It II 

Gun on Flat Substrate 




A1 Splat 

1 - 140 

2 

II II 

Gun on Flat Substrate 




Fe Splat 

1 - 100 

1 - 1000 

Ruhl (1967) 
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Table (3.2..L),- (contd.) 


System 

Sample Size 

h 

(cal/cm^ s°C) 

Source 

Gun on Flat Substrate 

Al-Cu Splat on Glass 
Piston & Anvil 

25 - 100 

0.956 

Scott (1974) 
Williams 

Al-Si Splat 

25 - 50 

24 - 215 

& Jones(1975) 

Gun on Flat Substrate 

A1 Splat 

100 

0.95 

Jones (1971) 

Gun Method Fe-Ni 

Glass Splat 

0.1 - 0.3 

24 - 240 

Davies(1978) 

Twin Roll Fe-Ni 

Glass Splat 

40 

2.4 - 24 

m ?i 

Chill Block MS 

Fe-Ni Glass Splat 
Chill Block MS 

20 

2.4 

it it 

Vincent 

Al-Si & Nimonic Splats 

20 - 50 

1.67 

et al (1980) 


Conventional Cast 
Al & Pb on Fe 

20000 

0.027 

Sully (1976) 

Continuous Cast 
Steel on Water 

50 000 

0.03 

Hills (1965) 

Cooled Cu 
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Table (3.2.1).- (contd.) 


System 

Sample Size 

h 

Source 


( y^m) 

(cal/cm^ s°C) 


Free Flight MS 
Metglass on Brine 

100 

0.04 

Kavesh(1976) 

Planar Flow MS 
Fe-B Glass 

20 - 40 

12 - 50 

Huang & 
Fiedler (1981) 

Atomization of 



Gill 

Undercooled A1 

50 

0.0678 

et al (1984) 

Melt Extraction 
Fe-Ni Wires 

25 - 1000 

0.1 

Robertson 
et al (1978) 

Direct Chill 
Horizontal 
Continuous Casting 
Al, Pb, Sn, and Zn 

20 000 

0.024 - 0.86 

Weckman & 

Niessen 

(1984) 

Atomization 

10 

2.4 

Cohen et al 
in Mehrabian 
et al (1980) 

Melt Spinning 

25 

2.4 

• i tr 


Self -Qu enc h ing 

10 

very large 

n it 

Rod Casting, Al 

50 000 

0.05 

Davies & 




Westby(1974) 
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The avoidance of the Newtonian cooling assumption can lead to 
additional , important results which cannot be obtained from the 
simpler models. Moreover, there is the hope that the fewer 
approximations that are introduced into the model the closer to 
reality its predictions will be. 

The opposite extreme to Newtonian cooling is to assume that the 
splat is in perfect thermal contact with the chill. This condition 
is known as ideal cooling. From the well known solutions to heat 
conduction problems under ideal cooling conditions (e.g. Carslaw 
and Jaeger (1959) ) and from Schwarz’s solution to the solidifica- 
tion problem (Sec(3A.2)), Jones derived approximate expressions 
for cooling and freezing rates under ideal cooling conditions. 
These are (see Jones(1982) ) : 

dT/dt = B/x 2 

and 

R = dx/dt = B'/x 

where x is the distahce in the slab ^from the chill and B and 
B' are functions of the relevant temperature intervals and of the 
material properties. 

It is perhaps not surprising that neither the assumption of 
Newtonian cooling nor the one of ideal cooling represented by 
Eqns(l) and (2) seem to be able to accurately represent observed 
behavior. This means that even though the thermal contact at the 


( 1 ) 

( 2 ) 
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splat— chill interface is far from perfect, the temperature gradients 
across the splat cannot be neglected. This, most frequently found 
cooling regime, is conveniently called intermediate cooling. 

In the sequel we present the details of the formulation and 
the solution of a model of a typical RSP system working in the 
intermediate cooling regime. In the description we will concentra- 
te on the PFMS system since the bulk of our calculations were 
performed for that configuration. However, summary comments on the 
modeling of other RSP systems are also presented. We hope our 
methods to be sufficiently general as to allow their application 
to any other RSP system. Thus, after a detailed description of 
the model of the PFMS system and of the results obtained from it, 
we discuss some points about the Twin Roll RS system, the Piston 
and Anvil system, the Melt Fragmentation processes and the systems 
based on Surface Heating and on Spray Deposition. 


3.3.1.- A Model of the Planar Flow Melt Spinning Process. 

The Melt Spinning (MS) process is one of the most commonly 
used methods of RST . The principle of the technique is very 
simple. A sample is melted inside a crucible and then a sudden 
pressure surge ia t applied to produce a thin liquid jet from a 
nozzle at the bottom of the crucible. This jet is in turn direc- 
ted towards the surface of a rapidly moving wheel. On impingement, 
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the liquid jet transforms into a small puddle , Finally, a thin solid 
ribbon is dragged from underneath the puddle by the moving wheel. 

Two main variants of the MS process exist, namely the Chill 
Block and the Planar Flow systems. The most significant difference 
between the two techniques is the detailed nature of the puddle 
formed at the point of impingement of the molten jet on the 
wheel. The arrangement used in the PF process restrains the puddle 
and promotes its stability. In Fig(l) we show schematic 
representations of the melt puddles formed, respectively , in the 
CBMS and in the PFMS process. It can be readily seen that in 
the latter the nozzle is brought into close proximity with the 
wheel . 

In the following pages we present our mathematical model of the 
PFMS process. We start by describing its formulation, then we 
discuss the solution methodology employed and conclude with a 
description of the results and some recommendations. 

a) Formulation of the Model. 

Because of its potential applications, the MS process has 
received a great deal of attention from the RS community. Most 
of the studies to date, however , have been concerned with the 0 
CB process whereas the interesting PF process has been left 
somewhat aside. Interestingly enough, the peculiar arrangement 
used in this latter system produces a fluid flow configuration 
very similar to those observed in Lubrication Technology 
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(a) 



Fig (3 . 3 . 1 . 1) Schematic comparison of melt puddle shapes for 
the (a) CBMS and (b) the PFMS systems. 
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systems. Since the two MS processes are closely related, we will 
start this section by reviewing previous modeling work in this area. 
Kavesh(1978) performed one of the first quantitative studies of 
the effects of changes in the process parameters on the properties 
of the products of CBMS . Using basic concepts from boundary 
layer theory he constructed representations of the fluid flow 
and the heat transfer phenomena taking place in the puddle. 

Based on the well known smallness of the Prandtl number values 
for metallic melts, he concluded that the mechanism responsible 
for the final ribbon thickness emerging from the puddle was the 
rate of heat transport out of the melt and into the chill. Kavesh 
also presented closed form expressions relating the geometry of 
the melt spun ribbons to important process parameters such as the 
melt flow rate and the wheel velocity. These relationships have 
been extensively used as a basis for many subsequent empirical 
studies of the process (e.g. Charter et al(1980)). 

Anthony and Cline (1978) also employed ideas from boundary layer 
theory in the vorticity-stream function formulation. They obtained 
closed form expressions demonstrating that the layer inside the 
puddle in which most of the temperature change takes place was 
many times larger than the shear layer due to wall induced 
vorticity. 

den Decker and Drevers(1980) used again boundary layer theory, 
this time combined with the equations of phase change kinetics by 
nucleation and growth and incorporating the temperature dependence 
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of the melt viscosity. They were able to compute both the final 
ribbon thickness and the crystal to glass ratio in the resulting 
ribbons as a function of the process parameters. In their conclu- 
sions, they agreed with the results of previous researchers in 
that the computed thermal boundary layers were much thicker than 
the corresponding momentum boundary layers. 

Katgerman (1980) modeled the flow in the puddle using the x— 
component of the momentum balance equation in a coordinate 
system fixed to the solidification interface. He also used 
approximate closed form expressions to represent the freezing 
and computed the thicknesses of the momentum and thermal boundary 
layers for several cooling conditions at the splat-wheel inter- 
face. He also concluded that the transfer of thermal energy 
played a predominant role in the determination of the final ribbon 
thickness . 

Vincent et al(1982) reviewed work on the modeling of MS and, 
using also concepts from boundary layer theory, concluded that 
momentum transport was the dominant mechanism controlling the 
final ribbon thickness. A similar conclusion was reached by 
Takeshita and Shingu(1983) who incorporated the temperature 
dependence of the melt viscosity in their calculations of flow 
and heat transfer using the equations of boundary layer theory. 

Three main points emerge as a result of the preceding review 
of the literature. First , although all investigators have used 
basically the same equations in one form or another, there is 
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no agreement as to the true mechanism controlling the final ribbon 
thickness. Momentum transport, heat transfer, and mixed mechanisms 
have all been proposed. Secondly , no distinction has been drawn 
between the CB and the PF processes , and the bulk of the 
simulations have been done for the first of these. We believe that 
the two processes , although very similar at first sight, possess 
however some features which warrant separate treatments. Finally 
we note that most of the previous models of the MS process have 
dealt with metallic glasses thus avoiding the problems associated 
with the release of the latent heat during solidification. 

Miyazawa and Szekely(1981 ,1979) have used a different approach 
to model other, somewhat related RSP systems. They have chosen to 
represent the flow phenomena taking place inside the splat during 
RS by a modified form of the Navier-Stokes equations in which 
the inertia forces are considered negligible compared with the 
viscous and pressure forces. The same approach has long been used 
by mechanical engineers to represent systems in which moving solid 
surfaces are separated by thin fluid layers and the resulting 
equations constitute the basis of the so called theory of 
lubrication. The use of lubrication theory to represent fluid 
flow in RSP systems was considered appropriate since small 
sample sizes in at least one spatial direction are one of the 
important features of these systems. In this sense they are 
analogous to the systems encountered in ball bearings. 

Encouraged by the results reported by Miyazawa and Szekely 


64 



for the piston and anvil and for the twin roll devices we decided 
to proceed further and investigate the use of lubrication theory 
for the representation of the flow behavior in the puddle of a 
PFMS system producing crystalline ribbons. 

The governing equations for fluid flow and heat transfer- 
solidification subject to the assumption of negligible inertia 
forces are (see Sec(3A.l)-(3A.3)) ; 


The equation of continuity 

div V = 0 

the equation of motion 



grad V) 


grad P 


the differential energy balance 


( 1 ), 


( 2 ), 


V ’ grad E 


div( K grad T ) 


(3), 


and the enthalpy-temperature relationship 


Ef + (T-I L )y»c p 
E f ((T - T s )/(T l - T g )) 
(T - T ) f C 


T > T 


T ^ T ^ T 

S L 


(4) 


T < T 
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Equations (1) — (3) may be simplified still further for the case 
of the schematic PFMS system shown in Fig(2) . So , by introducing 
the following simplifying assumptions, namely ; (i) planar flow 
conditions, (ii) negligible y-momentum compared with the x- 
momentum, (iii) velocity derivatives along the downstream 
coordinate negligible compared to derivatives across the puddle 
thickness, (iv)convectioh in the y-direction negligible compared 
to convection in the x— direction and viceversa for conduction, and 
(v)physical properties constant independent of temperature , 
Eqns(l)-(3) become, 

H 

Q/w = f V dy = (Q + Q,)/w (la) 

IX S -L 

O 


d 2 V /dy 2 = dP/dx 

X 


2 2 

V dE/dx = K d T/dy 

X 


Equations (la)-(3a) and (4) must now be solved subject to 
appropriate boundary conditions. The ones we have chosen , again 
based on the schematic of Fig (2) are; (i) partial slip or no- 
slip at the splat-wheel interface, (ii) melt flow rate constant 
and given, (iii) heat transfer at the splat-wheel interface 
specified by a heat transfer coefficient, (iv) heat flow through 
the top surface of puddle negligible compared to the heat flow 
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superalloy ribbons. The nozzle is at the top and 
the wheel at the bottom, moving from left to right. 
Courtesy of General Electric and NASA. 


through the wheel, and (v) at the melt-gas interface downstream, 
the shear stress is so small that can be neglected but the normal 
stress is determined by the capillary pressure given by Laplace's 
equation. 

We now briefly comment on the appropriateness of both the 
assumptions and the boundary conditions. Assumptions (i)-(iii) 
are considered adequate since the geometry of the system involves 
one characteristic spatial dimension which is much smaller than 
the other two. Assumption (iv) is plausible in view of the 
characteristic ability of molten metals to transmit heat more 
readily than momentum (small Prandtl number). Assumption (v) 
is also justified since we are interested in events happening 
inside the puddle and temperature extremes there are never 
greater thantwo or three hundred degrees. 

Boundary condition (i) is used mainly to remove the stress 
singularity that it is known to result from the use of the no- 
slip condition at the line of contact melt-gas-wheel upstream. 
Condition (ii) is appropriate since the melt flow rate is 
ultimately specified by the imposed pressure in the crucible. 
Condition (iii) is adopted mainly as a means of making up for 
our ignorance about the details of the complicated events taking 
place at the splat-wheel interface. However, the best available 
values of h have been used. Boundary condition (iv) is 
justified as can be readily checked by comparing the heat losses 
due to radiation from the top with the heat lost by contact with 
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the wheel. Finally, condition (v) has been adopted to be able 
to predict the exact location of the downstream meniscus and it 
is based on the ideas of capillary hydrodynamics as described, 
for example, by Levich(1962) . 

The complicating effects introduced by the consideration of 
capillary phenomena warrant further comment . As can be seen in 
Fig(2) , the top surface of the puddle detaches itself from the 
nozzle at some location along the downstream direction. After 
this, the top surface of the puddle becomes free in the sense 
that it is no more restrained by the solid nozzle lip. Moreover, 
this free boundary adopts a shape determined both by the capillary 
effect of the surface tension and by the solidification phenomena 
taking place underneath the puddle. We must note that, since the 
melt flow rate is constant, the amount of material passing through 
any given section at constant x is the same regardless of the 
particular location selected. 

The solution of the complete free surface problem is a complex 
matter. The addition of solidification effects would make the 
problem intractable if it were not for the introduction of suitable 
simplifying assumptions. Following Levich, we consider lubrication 
theory still valid after the detachment point and assume that the 
pressure acting across the puddle is given by Laplace's equation. 
After some manipulation (Sec(3A.3)), the following expression for 
the shape of the free surface is obtained , 


d 3 H L /dx 3 


(3 H 3 


xoyv) - 


Hi ) 


(5) 
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We have now completed the mathematical statement of the problem. 
What we have to do Is to solve Eqns(la)-(3a) , (4) and (5) simultaneous- 
ly subject to the boundary conditions mentioned above expressed in 
mathematical form. In the following section we describe our method 
of solution of this problem. 

b) Solution Procedure. 

We now describe the main features of the numerical algorithm 

we have used to solve the equations representing our model of the 

PFMS system. First, we note that for a given thickness of the 

solidified layer underneath the puddle, Eqns(la) and (2a) may be 

readily integrated to produce closed form expressions for V and 

x 

for P as function of the process parameters and of the material 

properties. However, in the region where the free surface forms 

after detachment from the nozzle, the thickness of the puddle is 

unknown in advance and must be calculated by solving Eqn(5) before 

we can proceed to integrate Eqns(la) and (2a) to get V and P . 

x 

Once V is known, its average across the puddle thickness 
can be calculated and then Eqns(3a) and (4) can be solved to give 
the cooling and freezing rates at every point in the puddle. Since 
no closed form analytical solutions exist for this solidification 
problem, we resort to numerical methods. The simplest explicit 
finite difference scheme has been selected to perform this 
calculation. The scheme can be proved to be stable and consistent 
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as long as the step lengths satisfy the stability condition, 

0 1 Ax/V x (Ay) 2 = 1/2 (6) 

Equation (5) must also be solved numerically. We have found 
that the transformation of Eqn(5) into an equivalent set of 
three ordinary differential equations of first order which are 
then solved using a simple Euler forward scheme to advance the 
solution in the downstream direction, produces an algorithm 
which is both stable and consistent. 

In summary, the complete set of steps we have used to solve 
the equations representing our model of the PFMS process is 
as follows: 

(i) Before the detachment point, Eqns(la) and (2a) are solved 

to obtain V and P . After detachment, however, Eqn(5) is 
x 

solved first, advancing one step in the downstream direction to 

find the location of the free surface and then Eqns(la) and (2a) 

are solved for V and P . In both cases , V is also 

x x 

computed . 

(ii) Equations (3a) and (4) are then solved to find T(x,y) and 
the location of the solidification interfase y 

s 

(iii) The puddle is swept in the downstream direction following 
the procedure indicated in (i) and (ii) until the free surface 
encounters the solidification interface. This is considered to 
be the end of the puddle. The point of intersection defines the 
final ribbon thickness. 
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We have constructed a FORTRAN computer program to perform 
the calculations indicated above. A typical run in the IBM 360 
computer of the IPS at MIT using 31 grid points in the y- 
direction and about 5000 in the x— direction, required 
approximately 60 s of CPU time. The main criterion for the 
acceptance of the results was the satisfaction of overall mass 
conservation and thus only the results of runs for which this 
was true down to less than 0.1% were accepted. The program 
contains many comments which should make it easier to understand 
it may readily be adapted into other computers having 
FORTRAN compilers. In the next section we describe and discuss 
the results that can be obtained with our program. The program 
listing itself, by the way, is included in Sec (5. 4) below. 


c) Description of Results and Discussion. 

We now present a summary description of the results that can 
be obtained with our model. The input data for the calculations 
have, for the most part, been kindly provided by G.E.-NASA. 
Additional data have been taken from the literature and, when 
v^l. u ®s were lacking, the best available estimates were used. 

The material under study was a Ni-base superalloy and the whole 
set of input data is shown in Table(l) below. 
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Table (3. 3. 1.1).- 
System. 


Input Data for the Calculations of the PFMS 


Nozzle Breadth 

0.064 

cm 

Crucible-Wheel Gap 

0.030 

cm 

Slot Width 

0.635 

cm 

Wheel Radius 

12.7 

cm 

Wheel Velocity 

1200 

rpm 

Ejection Pressure 

1.38 * 

10 3 g/cm 

Melt Flow Rate 

3.86 

3, 

cm /s 

Pouring Temperature 

1440 

°C 

Puddle Length 

0.29 

cm 

Ribbon Thickness 

0.0038 

cm 

Melt Density 

8.5 

g/cm 3 

Melt Viscosity 

0.046 

g/cm s 

Melt Surface Tension 

1778 

g/s 2 

Melt Specific Heat 

0.15 

cal/g °C 

Melt Thermal Conductivity 

0.0717 

cal/cm s 

Solidification Range 

1315 - 1335 °C 

Latent Heat of Fusion 

71.7 

cal / % 

Heat Transfer Coefficient 

1-2 

cal/cm^ 

Dynamic Contact Angle 

160 

o 

Initial Meniscus Curvature 

6.89 

cm 
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Figure (3) shows the computed temperature field corresponding 
to what we have called the typical data set (see Table (1)). It 
should be noted that in Fig (3), as well as in Figs (2) and (4), the 
vertical scale has been enlarged considerably for clarity of 
representation. In reality, for the horizontal dimension shown, 
the vertical scale is about 1/10 to 1/20 of the size shown in 
these figures. This shoul give an idea of the size scales in- 
volved in the problem. Moreover, this should make clear that the 
isotherms in reality lie almost parallel to the surface of the 
substrate, and that most of the temperature change takes place 
inside a relatively thin layer close to the bottom surface of the 
puddle. 

Other noticeable features in Fig (3) are as follows. First, as 
can be expected from the assumption of imperfect thermal contact 
at the splat-wheel interface, freezing does not start immediately 
on impingement at the point where the gas , the melt, and the 
wheel meet forming a line of contact. Moreover, when instead of 
a heat transfer coefficient we use a closed form expression 
representing ideal thermal contact (i.e. h ©o ) , solidification 
always started at the contact line. The resulting ribbon 
thicknesses , however, proved to be in all cases physically 
unrealistic since the solidified layer grew apparently too fast. 
Consistently throughout this work, the physically more satisfac- 
tory results were obtained when the assumption of non-ideal 
thermal contact at the splat-wheel interface was used. 
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Fig(3 .3 .1 .3) Computed isotherm map inside the melt puddle. The 
aolidifucation interface and the computed meniscus are also shown. 


Also seen in Fig (3) is the inflexion point appearing on the 
isotherms as they approach the free surface. This change in 
curvature , corresponding to an increased cooling rate, coincides 
with the formation of the free surface after detachment. Apparent- 
ly, the continuously diminishing puddle thickness observed after 
detachment has an influence on the thermal behavior in the form 
of increasing cooling and freezing rates. This observation leads 
to predictions about the variation of microstructure across the 
ribbon thickness which are in reasonably good agreement with 
measurements as described in more detail below. 

Figure (4) illustrates the computed streamline pattern 
corresponding to the same typical run. We recall that when the 
values of the stream function on every pair of neighbouring 
streamlines differ by the same amount from one pair to the other, 
it is possible to perceive the variations in the magnitude and 
the direction of the velocity over the domain at a glance since 
the same mass is flowing between any pair of such lines. The 
existence of two characteristically different flow domains inside 
the puddle is readily noted. The closely spaced streamlines 
close to the moving wheel indicate the existence of large 
splat velocities there. The circulating streamlines lying on top, 
on the other hand, show that there is a large region of slowly 
recirculating flow centered approximately in the middle of the 
puddle at about the point of detachment. 
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3.3.1.2_, - Photographs of the longitudinal cross 
section, (a), and the top surface, (b)-(d), of 
planar flow melt spun ribbons obtained by intro- 
ducing refractory powder together with the melt. 
The powder agglomerated and formed lumps inside 
the puddle. These lumps came out from the puddle 
at regular intervals of time. From Zielinski, P.G. 
and D.G. Ast(1983). 


/ ? 



These regions of , respectively, recirculating and forward 
moving flow are separated by the streamline marked >■- 1.0. 

This streamline , in turn, intersects the free surface at some 
location downstream forming what is called a stagnation point 
characterized by the velocity having the value of zero. It can 
also be seen in this figure that, inmediately after detachment, 
the layer of fluid lying above the partially solidified ribbon 
is dragged forward by the latter. However because of the melt 
fluidity, the liquid layer continues to thin down until the 
final ribbon thickness is reached. A final point worth noticing 
is the intense motion generated near the top surface of the 
puddle upon detachment resulting from the disappearance of the 
constraining effect of the nozzle wall. 

The characteristic fluid flow field computed by our model 
certainly represents an alternative picture of the system when 
compared with the results of all previous mathematical studies 
of the flow in the melt puddle. Although the upward curvature 
of the streamlines near the wheel had been predicted before, and 
explained by the increasingly larger portion of the total flow 
rate carried by the partially solidified ribbon underneath moving 
downstream, no one had found a region of recirculating flow. 
Despite the difficulties involved in resolving this point we 
note that there seems to be some empirical evidence which 
strongly suggests the existence of such recirculatory flow in 
the puddle. Mention can be made, for example, of the experiments 


80 


performed by Zielinski and Ast(1983) where finelly divided powder 
was introduced together with the melt. They suggest that the 
regularly spaced lumps observed on the top surfaces of ribbons 
produced by the PFMS process resulted from the agglomeration 
of individual powder particles during the intense recirculating 
motion taking place inside the puddle. In this sense, our results 
provide the first quantitative evidence for the existence of 
recirculatory motion inside the puddle of the PFMS system. 

The photographs showing Zielinski and Ast f s results have been 
included here in Plate (3. 3. 1.2). 

We note , finally , that the main features of the flow 
phenomena taking place inside the PFMS puddle have also been 
observed in other , somewhat related flow systems, namely in 
the study of lubrication with cavitation (Savage(1977) ) and in 
the analysis of coating flows (Kistler and Scriven(1984) ) . 

Figure (5) is a plot of the average cross-sectional cooling 
rate calculated according to the formula 

_ H 

T = (1/H) f V x (dT/dx) dy (7) 

0 

against the downstream coordinate. One can readily see that the 
cooling rate approaches 10^ °C/s inside the puddle before the 
start of solidification. However, the cooling rate peaks and then 
decreases, first gradually and faster later as freezing sets in 
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according to Eqn(7). 



and the latent heat of fusion is released. Moreover , once the 
solidification is complete, and in agreement with the results of 
calculations for the Twin Roll system, the cooling rate rises 
again, this time reachinga value above 10 6 °C/s, mainly because 

of the smallness of the section being cooled. Most likely, since 
the splat is not firmly attached to the wheel, once the ribbon 
leaves the puddle the heat transfer coefficient should decrease 
and these high cooling rates may only be rarely observed. 

It may be worth mentioning that a plot similar to Fig (5) 
can be constructed using the much simpler lumped parameter models 
described in Sec (3. 2). So one may justifiedly ask why is it 
necessary to use the more complicated models which take into 
account the temperature gradients across the splat in order to 
represent mathematically the PFMS process ? .This can be 
answered by recalling the main purpose of our work which is 
the establishment of the relationships between the process 
parameters and the structure-properties of the resulting product. 
It should be understood that a model which neglects the existence 
of temperature differences across the splat will be unable to 
predict any microstructure variations across such thickness. 

To point at the significance of the temperature differences 
existing across the thickness of the splat, in Fig (6) we present 
the computed temperatures on the top and bottom surfaces of the 
puddle/splat formed in the PFMS process, as a function of the 
downstream coordinate. As expected, one can readily see the 
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Where (dT/dy) is the temperature gradient in the (approximately) 

normal direction to the solidification interface, (dy /dx) is 

s 

the rate of thickening of the solidified layer along the down- 
stream coordinate and is the x- component of velocity at 

the interface. From Eqn(8) , can be calculated as a function 

of the downstream location. To compare our calculations with 
available experimental data we have used the fraction solidified 
(relative to the final ribbon thickness) instead of the x 
coordinate. The result of this calculation is shown in Fig(7). 

The very large values of the cooling rate prevalent during the 
first 1/10 fraction solidified are readily apparent. Also 
clearly seen is the very sharp decrease from the large initial 
values. A somewhat unexpected feature is the hump appearing 
arounf f = 0.3 . This hump coincides with - and it may be the 

result of - the detachment of the top surface of the puddle from 
the nozzle. Finally, as the end of solidification is approached, 
the effective cooling rate goes under 10^ °C/s . 

Considerable variations in the size of microstructural features 
have long been observed in melt spun ribbons . For the alloys of 
this study, measurements have been made to quantify this feature. 

A photomicrograph of a longitudinal section of a ribbon has been 
included here in Plate (3. 3. 1.3). We then decided to use the 
result of our thermal calculations to attempt a prediction of 
these microstructural variations. We started by looking for a 
suitable dendrite arm spacing-cooling rate correlation. We were 
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very rapid decrease of the temperature of the bottom surface of the 
puddle following impingement. Somewhat unexpected, however, is the 
much slower decrease of the temperature of the top surface. Moreover, 
when the cooling rate on the bottom surface slows down as a 
consequence of the onset of solidification, the rate on the top 
surface increases markedly. This increase in cooling rate coincides 
and it may be related to the detachment from the nozzle. Finally, 
similarly to what happens to points along the bottom surface, the 
cooling rate slows down as the free surface approaches the solid- 
liquid interface. It can be seen that the cooling rate during 
solidification is somewhat smaller for the top surface than for 
the bottom surface of the ribbon. This is to be expected from 
the relative location of the chill with respect to these two 
surfaces .Finally , mention should be made also of the considerable 
temperature gradients existing in the puddle. We can expect 
these appreciable temperature differences (reaching even thousands 
of °C/cm )to have a measurable effect on the nature of the 
resulting microstructure. 

To further investigate this point we decided to compute the 
"effective" cooling rate experienced by points just ahead of the 
solidification interface. The effective cooling rate was 
computed from the following formula, 

T e = (dT/dy) (dy g /dx) V x (8) . 


85 




not able to find the exact correlation valid for the alloys of 
our study so we decided to use a relationship proposed for the 
superalloy Inconel 718; this is (see Sec(2A.l)), 

. -0.34 

A = 34 ( T ) (9) 

s e 

The result of combining Fig(7) with Eqn(9) is shown in Fig(8). 

In the same figure also appear the results of actual metallographic 
measurements performed by Ms Segal (1983) at MIT. The plot just 
shows secondary dendrite arm/cell spacings as a function of the 
fractional distance from the wheel surface (referred to the 
final ribbon thickness) . The dashed line describing our results 
overestimates the spacing for the first third of the ribbon but 
corrects itself afterwards and remains in good agreement with the 
observed values during the last 2/3 of solidification. At this 
point we can only speculate on the reasons for the observed 
discrepancy during the initial stages of solidification. Several 
possible explanations can be offered. Among these we may mention: 

(i) The difficulty of accurately computing the value of the 
effective cooling rate at the beginning of solidification. The 
size of our computational grid being the ultimate limitation. 

This could produce a smaller rate of decrease of the effective 
cooling rate during the initial stages of freezing thus making 
the computed dendrite size curve closer to the measured values, 

(ii) the existence of phenomena unaccounted for by our model, as 
may be the case of undercooling in the melt layer next to the 
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3. 3. 1,^3.- Longitudinal section of melt spun 
ribbon of Inconel 718. The microstructure is 
similar to the one found in conventional castings, 
with the cell (dendrite) size increasing from 
bottom (wheel side) to top (free surface) . From 
Warrington, D.H. et al., in Masumoto, T. and K. 
Suzuki (eds) , (1982). 
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chill; the microstructure resulting from the solidification of 
such an undercooled layer can be expected to differ considerably 
from the ones characteristic of smaller cooling rates, (iii) it 
can also be the case that we just simply used the wrong dendrite 
size-cooling rate correlation. However, despite the apparent 
disagreement it is indeed remarkable that our model is capable of 
predicting microstructural sizes and cross sectional variations 
of microstructure of the correct order of magnitude, despite all 
the assumptions and uncertainties involved in our calculations. 

One important feature of our model is the coupling of the 
solidification phenomena occurring next to the wheel surface to 
the capillary processes taking place on the free surface. Thus, 
one of the products of our calculations is always the precise 
location of the melt-gas interface in the downstream side. In 
Fig (9) we show a plot of the location of such meniscus as a 
function of the roll velocity. One can readily note that the 
three computed free surface shapes are all very similar and indeed 
almost identical throughout the first half of the puddle. We 
note that the same detachment point was used in these calculations. 
However, as the free surface approaches the solidification 
interface, the effect of freezing on the shape of the meniscus 
becomes more significant. Since the total flow rate is constant, 
for any given downstream location, the sum of the partial flow 
rates carried by , respectively, the solid ribbon and the liquid 
film on top of it, must be equal to the total flow rate. So, 


91 






although the solidification rate is not altered significantly for 
changes in the wheel velocity within the range considered here, the 
total mass flow rate is an independent process parameter controlled 
mainly through the ejection pressure. Thus , for any given down- 
stream location, the partially solidified ribbon carries with it 
a larger fraction of the total flow rate when the roll velocity is 
large. This , in turn, has the effect of pulling the free surface 
down producing a thinner liquid film. 

The behavior described above contrasts with the one observed 
in the somewhat analogous system obtained when a solid object is 
withdrawn from a liquid bath. In this case a thin layer of fluid 
adheres to the surface of the object. The thickness of this film 
decreases with increasing distance from the exit point from the 
bath until it reaches a limiting value under steady state con- 
ditions. The final liquid film thickness is calculated to be 
proportional to some fractional positive power of the withdrawal 
speed. However, when withdrawing solids from liquid baths, the 
supply of liquid is practically unlimited and the solid body is 
able to carry as much fluid as it can. On the other hand, during 
melt spinning, the melt flow rate is given by the ejection 
pressure and the final ribbon thickness represents a compromise 
between the imposed flow rate and the chilling effect of the 
wheel , 

An important feature of the PFMS system is the widespread 
use of slot shaped nozzles. These rectangular nozzles allow for 
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the production of ribbons or sheets with a width essentially equal 
to the width of the slot. For comparison, in the CBMS process. 


circular nozzles have been used most of the time. The cylindrical 

jet produced during CBMS , on impingement, spreads laterally 

on the moving wheel producing ribbon with a width which is several 

times larger than the diameter of the initial jet. Since the 

resulting spreading phenomenon is difficult to control, variations 

in width along the ribbon length are commonly observed. On the 

other hand, during PFMS , the ribbon width is determined by the 

width of the nozzle and since both the melt flow rate and the 

exit velocity of the ribbon ( = V ) are fixed, an overall 

r x 

mass balance leads to 


H = Q/(wV ) (10) 

r r x 

The behavior represented by Eqn(10), remarkable because of its 
simplicity, is in sharp contrast with the more complicated 
relationship involving the same quantities found to be valid in 
the case of CBMS (see e.g. Kavesh) . Needless to say, the more 
complicated relationship in the case of CBMS arises because 
of the additional factors involved with the lateral spreading of 
the molten jet on impingement on the wheel, a phenomenon which 
is characteristically absent during PFMS with slot shaped 
nozzles. From Eqn(10) we see that we should expect the final 
ribbon thickness to be inversely proportional to the wheel 
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velocity. We decided to run our program several times maintaining 
everything constant except the wheel velocity (and, of course, the 
melt flow rate) . The result of such an experiment is shown in 
Fig(10). As expected, the results of our calculations fall very 
well along a straight line of negative slope. So, our results 
certainly satisfy the overall mass conservation requirement. This, 
however, is to be expected since the overall mass balance was used 
throughout as a criterion for the correctness of the results. 

It is indeed unfortunate that very few experimental data of this 
kind, produced under carefully controlled conditions, are 
available. The only experimental data we got , for the system 
of this study, are included in Fig(10), However, caution must be 
used when comparing measurements and calculations in this plot 
since all the experimental data points were obtained under 
entirely different process conditions, in particular , under 
different melt flow rates. Nevertheless, it is encouraging to see 
that our model predicts expected trends accurately. 

Furthermore, we should note that, once the wheel velocity and 
the melt flow rate are specified, the computed line shown in 
Fig (10) is related to one and only one value of the heat transfer 
coefficient at the splat-wheel interface. We could then expect 
empirical curves analogous to the theoretical one shown, be used, 
in conjunction with a model like ours, for the determination of 
h . 
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the final ribbon thickness vs wheel velocity relationship. 



From Eqn(10) we should also expect the final ribbon thickness 
to be directly proportional to the melt flow rate. We decided to 
perform several runs of our program in which the only parametric 
variation was the flow rate. Perhaps not surprisingly, we found 
that the overall mass balance could not be satisfied unless the 
length of the puddle was adjusted accordingly. This change in 
the size of the puddle with the flow rate is not unique to our 
calculations but has also been consistently found in the 
laboratory (see e.g. Huang and Fiedler (1981) ) . The results of 
our calculations of the influence of the melt flow rate on the 
final ribbon thickness and on the puddle length are shown in Fig 
(11) . The results suggest that thicker ribbons may be produced 
simply by increasing the melt flow rate. This conclusion may be 
deceiving since larger flow rates produce bigger puddles which 
in turn are increasingly unstable because of surface tension 
effects. Moreover, as described before, thicker tapes would be 
subjected to larger temperature differences across their thick- 
nesses unless something is done to attenuate the chilling effect 
of the wheel. The resulting thermal history will invariably lead 
to more significant variations of the microstructure of the strip 
across its thickness. Furthermore, it is clear that larger flow 
rates, because of the attending increased energy content, will 
require the substrate to be able to absorb more heat to be able 
to induce rapid solidification on the splat. The expert’s 
consensus seems to be that melt spinning wheels are operating at 
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Fig (3.3.1.11) Predicted effects of the melt flow rate on 
the final ribbon thickness, and on the size of the puddle 



or very close to their maximum heat extraction capacities. 

The heat extraction ability of the melt spinning wheel may be 
well described, in the context of this work, by the corresponding 
value of the heat transfer coefficient at the splat-wheel 
interface under given operating conditions. It must be clear that 
the value of such coefficient does not only depend on the physical 
properties of the wheel and the splat but it is intimately related 
to the actual values of the various process parameters. To further 
investigate this point we performed some calculations changing the 
value of h and adjusting the melt flow rate accordingly in 
such a way that the overall mass balance was always satisfied. As 
expected, the heat transfer coefficient appeared to have a strong 
influence on the final ribbon thickness achieved. The actual result 
of this calculation is shown in Fig (12). The point mentioned before 
can be seen here more clearly. Much more efficient heat absorbing 
substrates are required in order to produce melt spun ribbons 
using large melt flow rates. This constraint is so strong that in 
reality it is practically impossible to induce high cooling rates 
(of the order of , say 10“* °C/s) in ribbons thicker than, say 
one hundred microns. 

The superheat imposed on the melt prior to spinning is one 
easily controlled parameter. Figure (13) shows the computed final 
ribbon thickness as a function of the superheat of the melt. As 
expected on simple physical grounds, the final thickness decreases 
as the superheat is increased. This effect is particularly intense 
at high superheats. However, this results must be Interpreted with 
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Fig (3.3.1.12) Predicted effect of the heat transfer 
coefficient at the splat-wheel interface on the final ribbon 
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caution since the property changes produced by the higher pouring 
temperatures can modify the wetting behavior of the melt and lead 
to the apparently contradictory findings of some researchers 
(see e.g. Scott (1974) ) . 

d) Conclusions 

1) We have constructed and programmed a mathematical model 
of the PFMS process which is capable of predicting the dynamic 
behavior of the system, particularly regarding the relationship 
between process parameters and ribbon microstructure. Comparison 
of the computed results with the limited empirical data available, 
together with the extensive checks contained in our program, point 
towards the correctness of our approach. 

2) From our simulations we conclude that the most important 
process parameters of the PFMS technique are: 

(i) The melt flow rate, 

(ii) the wheel velocity, 

(iii) the heat transfer coefficient, 

(iv) the geometry of the system, and 

(v) the physical properties of the materials involved. 
Interestingly enough, a similar set of relevant process parameters 
was arrived at empirically by the NASA group (Jech et al(1984)). 

Only a restricted set of values of these process parameters 
allows for the steady production of uniform ribbons. The restrict 
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tion to particular sets of values of the process variables, although 
significant, seems less stringent than the one observed in the case 
of the Twin Roll system. There seems to be an advantage in using the 
single roll device for RSP in view of its increased flexibility 
of operation and of the possibility for continuous processing. It 
would seem that the inevitable deterioration of both wheel and 
nozzle due to thermomechamical stresses and chemical effects is 
one of the main obstacles to the continuous production of strip. 

3) We have demonstrated how the basic principles of lubri- 
cation theory, capillary hydrodynamics and solidification heat 
transfer can be combined to produce a meaningful picture of a 
typical RSP system. The insight gained from this approach 
suggests the possibility of useful information to be obtained 
by extending our methods to other RSP systems. 

4) We believe that a computer program such as the one 
described here may well constitute the core of an on-line control 
mechanism for the automation of the PFMS process. Ideally, 
however, the program should be fine-tuned by comparing its 
predictions with the results of a carefully designed and con- 
trolled set of experiments in which not only the ribbon is 
characterized in terms of its microstructure and properties but 
all the relevant process parameters are also measured or 
controlled as accurately as possible. 

5) The FORTRAN listing of the program used to perform the 
calculations described in this section is included in Sec(5.4). 
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3.3.2.- 


Some Comments on the Modeling of Other RSP Systems. 


In the next few pages we present a bird's eye view of some 
important aspects of the mathematical modeling of some selected 
RSP systems. Examples from all groups presented in Table (3. 1.1) 
will be discussed. However, in all cases the discussion will be 
brief and limited to a few important points. We start with the 
Twin Roll and the Piston and Anvil RS processes since they are 
both good examples of splatting methods and are also somewhat 
related to the PFMS system described in the previous section. 

We then move on to the Melt Fragmentation techniques and summary 
comments are included about the various atomization processes. 
Finally, mention will be made of the modeling of RS Laser 
Processing and of Spray Forming techniques. 

a) Modeling the Twin Roll RS Process. 

The most comprehensive model of the TRRS process to date has 
been presented by Miyazawa and Szekely(1981) . Starting from the 
momentum balance equations in the lubrication approximation, they 
solved the differential energy balance equation with due account 
taken of the latent heat of solidification. Besides the thermal 
calculations, closed form expressions were obtained for the 
velocity field of the material in the roll gap. The phenomenon of 
solid deformation was accounted for by considering the solidified 
strip to be a creeping solid changing shape according to the 
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Norton-Bailey law of secondary creep. 

For the solution of the thermal problem, an implicit front- 
tracking method was used after a coordinate transformation of 
the original equations. The pressure profile in the roll gap was 
obtained by integrating the expression for the pressure gradient 
along the casting direction using a Runge-Kutta method. We have 
repeated the pressure calculations performed by Miyazawa and 
Szekely but this time including the material parameters describing 
the strain rate in the solidified ribbon as given by Frost and 
Ashby(1982) . Very large peak values of the pressure were encoun- 
tered for the conditions used (Fig(l)). The listing of the 
program used to perform this calculation is included in Sec(5.5). 

Among the most important conclusions of that study are the 
following: 

(i) Only a very narrow range of values of the process parameters 
allows for the steady operation of the process. Experimental 
precautions are required in order to obtain satisfactory perform- 
ance. This has been verified in practice (Murty and Adler (1982) ) . 

(ii) The main parameters of the technique are; the roll gap, the 
roll velocity, the feed rate and the material properties. The 
final ribbon thickness decreases with increasing the roll speed 
and by decreasing the roll gap or the feed rate. 

(iii) The roll separating force increases with decreasing the 
roll speed because of the increased reduction ratio. Moreover, 
slip between roll and splat accounts for the experimentally 
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observed differences in velocity of wheel and ribbon. 

Miyazawa et al(1983), have continued their study of the TRRS 
process. These papers should be consulted for further details. 


b) Modeling the Piston and Anvil System. 

Bletry(1973) presented the first detailed calculations of 
temperatures and cooling rates for the piston and anvil system. 
The next important step was performed by Miyazawa and Szekely 
(1979) who incorporated the spreading and squeezing phenomena 
into the thermal calculations. They used the momentum equations 
in the squeeze flow approximation to obtain expressions for the 
velocity field. These were used , in turn, together with an 
explicit finite difference form of the differential energy 
balance equation, to calculate the cooling and freezing rates 
for the process. The most important process parameters of the 
technique were found to be; the piston speed, the sample size, 
and the pouring temperature. The main variables influenced by 
the values of these parameters are the final splat thickness 
and the attainable cooling rates. 


c) Modeling the Melt Fragmentation Techniques. 

The two main groups are , as indicated in Table(3 . 1 , 1) , the 
atomization processes and the impact disintegration processes. 
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Mathematical modeling techniques have been used only very little 
in the study of these processes. The explanation for this may 
be the highly complex nature of the phenomena which take place. 

We limit ourselves here to mentioning the main variables 
involved and the literature sources where more information can 
be found. 

In centrifugal atomization, a spinning dish, cup or electrode 
produces fine particulates by first thinning of the melt and 
subsequent breaking up of the thinned layer. The main process 
parameters of this process have been found to be (see e.g. 

Hodkin et al(1973), Schmitt (1979) and Champagne and Angers (1984) ) ; 
the dish diameter and speed, the melt density, surface tension, and 
viscosity, the feed rate and the pouring temperature. Using basic 
concepts from fluid mechanics, relationships have been derived 
giving, for example , the fluid film thickness at the rim of the 
rotating dish and the resulting mean droplet diameter. Interestingly 
enough, many similarities with the well known centrifugal atomiza- 
tion process widely used in chemical engineering technology, have 
been found. There is one additional complicating factor in this 
case, however, since the fluid not only has to be fragmented but 
the resulting particles must be cooled and solidified. 

Gas and water atomization processes have long been used by the 
powder metallurgy industry. The introduction of ultrasonic gas 
atomization , with its resulting smaller particle size, has carried 
the fluid atomization processes into the realm of RST . The main 
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process parameters in this case have been recognized to be; the 
molten metal properties and flow geometry and the gas (water) jet 
flow rate and geometry. Fluid dynamic phenomena are strongly 
involved in the thinning and breaking up of the liquid metal 
jet during fluid atomization. The heat transfer processes of 
cooling and solidification which take place simultaneously add 
considerable complexity to the system. Tallmadge, for example, 
(1978) has performed extensive studies of metal atomization 
processes from a chemical engineering point of view, whereas 
Beddow(1978) and Lawley(1977) have presented comprehensive 
reviews of the subject. Grant (1983) also presents an overview 
of atomization processes but from the perspective of RST . 

The calculation of the cooling and freezing rates of the 
undercooled metal droplets produced by metal atomization 
processes has been performed most recently by Gill et al(1984). 

In the case of the vibrating electrode process (Ruthardt and 
Lierke(1981)) the melt is fragmented by the break up of the 
capillary waves induced on the thin melt layer by the vibrating 
electrode, In this case, the theory of capillary waves has 
been used to estimate the mean droplet size as a function of the 
melt surface tension, density, and viscosity and of the 
vibration frequency of the electrode. 

Schmitt (1979) has described the basic physics of impact dis- 
integration processes. He found that the final particle size in 
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this system depended on ; the nozzle diameter, the impactor radius 
and size, the distance from the nozzle to the impact point, the 
impactor rotation velocity, and the material properties. Lynch 
(1982) has used the basic principle of impact disintegration 
processes to produce particulates from a continuous molten metal 
jet. 

d) The Modeling of RS Laser and Electron Beam Processing and 
of Liquid Dynamic Compaction (Spray Forming) . 

The scanning of the surfaces of bulk materials with high 
intensity laser or electron beams can produce rapid melting and 
solidification of thin surface layers. In these processes the most 
important parameters are; the wavelength of the radiation, the 
incident power density, the interaction time , the detailed 
nature of the surface and the material properties. 

Breinan and Rear (1983) have modeled RS laser processing by 
using the one dimensional heat conduction equation with a source 
term. In this way they were able to avoid the complicating 
factors introduced by the change of phase and the fluid motion 
inside the melt layer. They claim very good agreement between the 
results of their calculations and those obtained from a more 
sophisticated finite element analysis which took the latent heat 
into account. Moreover, they felt that the accuracy of their 
model was good enough for comparison with their own measurements. 
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Chan et al(1983) have presented a more comprehensive model of 
IS Laser Processing in which they take into account the fluid 
flow phenomena taking place in the molten layer. The heat 
transfer and fluid flow problems are solved simultaneously 
and predictions are made about cooling and freezing rates as 
well as of puddle morphology. Mawella (1984) constructed a 
thermal model for the Electron Beam Processing of solids 
based on the theory of surface sources of heat, popular in 
welding calculations . He claims good agreement between his 
calculations and the results of his own experiments. 

Liquid Dynamic Compaction (also Spray Forming or Spray 
Deposition) is the name given to a group of RS processes 
in which relatively thick sections of rapidly solidified material 
are prepared by the continued showering of a substrate with a 
spray of molten metal droplets (which in turn may be produced by 
any of the various atomization processes). The structures and 
properties of the resulting products are functions of both the 
spray characteristics and the properties of the substrate. Singer 
and Evans (1983) have described the fundamentals of a statistical 
model for the representation of the spray. They have also suggested 
how to use the heat flow equations for the calculation of the 
specific conditions leading to the steady operation of the process. 
They have found the predictions of their model to be in good 
agreement with measurements performed in their own laboratory. 

Apelian et al(1983), and El-Kaddah et al(1983) have performed 
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detailed heat transfer calculations for the closely related Plasma 
Deposition processes. They found that the events taking place in 
the plasma are as important or more that those happening once the 
molten particles impact the substrate. Their papers should be 
consulted for details. 

In summary , although a large amount of work has been done 
on the mathematical modeling of RSP systems, there seem to be 
ample opportunities for additional research . Many systems have 
been studied only partially and nothing at all is known about a 
few others. Because of the relatively novel nature of RST , 
mathematical modeling can contribute to substitute costly trial 
and error procedures by more rational approaches to process 
development and control . 


112 


Chapter 4 


CONCLUSIONS AND SUGGESTIONS FOR FURTHER WORK 


Since most of the conclusions pertaining to the modeling of the 
PFMS process were presented before (see Sec (3. 3.1)), in this 
chapter we only add a few more points not explicitly mentioned 
there. On the other hand, we will include some conclusions of a 
more general nature about the potential of the mathematical 
approach to help in the understanding of the complex features 
of RSP systems. Specific suggestions are made regarding the 
directions in which mathematical modeling can be applied to the 
study of RST . In the sequel we simply present our points 
without following any particular order. 

a) One important limitation of the PFMS process is the result 
of two conflicting requirements. On the one hand, to obtain RS 
effects a small sample size is needed in at least one spatial 
direction. On the other, the production of thicker strip requires 
larger melt puddles. Since surface forces are what holds the 
puddle together during the production of thin ribbons, as mentioned 
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in Sec(3 .3.1) above, it is clear that the question of the stability 
of the melt puddle is critical to the success of the process. It 
may be seen that the production of relatively thick tapes is faced 
with big problems. It might be helpful to artificially constrain 
the puddle in some way instead of letting the surface forces 
do the job alone. However, according to our results, the thicker 
the produced tapes, the more significant the microstructure 
variations across the splat thickness will be. It is an intrinsic 
feature of the heat transfer-solidification processes taking 
place in this system that causes the highest cooling rates to 
be limited to that portion of the splat which is closer to the 
moving chill . 

b) It is well known that the unavoidable roughness always 
present on the wheel surface during PFMS is one important 
reason for the considerable variations in surface quality and 
microstructure of melt spun ribbons. The uneven wheel surface 
produces localized regions of relatively good thermal contact, 
which are separated by areas where the contact is poor (lift-off 
areas) . Ways must be devised to deal with this problem if one 
wishes to optimize the process. One possible alternative frequently 
suggested consists in homogenizing the cooling power of the wheel 
surface by eliminating the areas of good thermal contact. One way 
of doing this is by coating the surface of the wheel by a thin layer 
of a heat insulating material, (e.g. glass). This layer should be 
sufficiently thick to thermally homogenize the surface but thin 
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enough as to not to impair the heat extraction ability of the wheel. 


c) More attention should be paid to the spreading and wetting 
phenomena taking place during PFMS. In particular, the contact 
line formed in the rear of the puddle at the point of impingement 
should be looked at more closely. This contact line is the locus 
of many important processes involving fluid flow and heat 
transfer and the understanding of these may prove to be vital for 
the successful implementation of the technique. Both experimental 
and theoretical work is required in this area. 

d) The additional effects introduced by the presence of inertia 
forces (which were neglected in our calculations) should be more 
fully explored. Although because of the large surface tensions 
characteristic of metallic melts, this effects are likely to be 
small, more work is needed to verify this expectation. 

e) An effort should be made to perform carefully planned 
experiments aiming at the verification of the existing mathematical 
models of RSP systems. Work can proceed along the lines of 

the research performed on the CBMS system or following the 
pattern set by Miyazawa et al.. It is indeed unfortunate that 
so little information of this type is available for the PFMS 
process. In any case, the aim of experimental programs should be 
the establishment of the relationships linking the process 
parameters to the structure and properties of the resulting 
products. 
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f) The entire area of melt fragmentation processes remains a 
challenging ground for modeling work. The powder metallurgy 
industry has motivated much of the existing work on metal atomi- 
zation. However very little of this research has been of the 
theoretical kind. Many aspects of the various melt fragmentation 
techniques still remain obscure with empirical correlations being 
the only quantitative means of studying the processes. There 
are considerable opportunities for useful contributions from 
mathematical modeling in this whole area. The abundance of complex 
hydrodynamic phenomena coupled to heat transfer-solidification 
processes always present in these systems should attract the 
attention of the mathematically oriented engineer as well as that 
of the fluid dynamicist or applied mathematician. Work in this 
area could start from the important contributions by Hinze(1955), 
Dombrowski and Johns(1963), and , particularly, the one by 
Bradley (1973) . A good summary of the theory of fluid jets and 
their stability can also be found in Anno(1977). 

g) A somewhat related set of problems is found in the area of 
liquid dynamic compaction (spray forming). Here, the detailed 
structure and properties of the deposited product depend on the 
fluid dynamics of the swarm of droplets of varying sizes, on the 
interaction of such droplets with the underlying substrate and on 
the heat transfer phenomena taking place both during flight and 
after impingement of the droplets on the deposit. 
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Needless to say, obtaining deposits with low porosity and desirable 
properties should be achieved more easily by using a rational 
approach in which the use of empirical trial and error is minimized. 

h) The use of high intensity energy beams to alter the surface 
the bulk samples by rapid heating and cooling is now well establish- 
ed at an industrial level. Heat treating, welding, surface melting 
and freezing and surface alloying have all been demonstrated. 
However, despite the successes, there are still areas where 
knowledge is scanty. 

In all surface heating systems one invariably finds thin layers 
of material subjected to very large temperature gradients. In 
particular, during RS laser or electron beam processing, thin 
molten layers form on the surface of the object which solidify 
once the heat source is removed. The behavior of these thin layers 
of fluid is governed by the same basic principles of continuum 
mechanics. Because of the small scale, capillary phenomena play an 
important role in these systems. We believe that increased under- 
standing of the detailed behavior of these systems can be obtained 
from a more extensive application of the equations of fluid 
dynamics and heat transfer with change of phase. Even more know- 
ledge can be gained if these thermal calculations are coupled with 
crystal growth kinetics for the prediction of the microstructure 
of the surface layers resulting from processing. 
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i) The intrinsic limitation of RSP to small sample sizes invary- 
ably results in systems with a disproportionately large surface 
area to volume ratio. This makes imperative a more detailed study 

of the interfacial phenomena taking place at interphase interfaces. 
The peculiar behavior of lines of contact, which are almost always 
present in RSP systems certainly deserves more careful consid- 
eration. Very often, these contact lines are the locus of complex 
heat transfer and fluid flow phenomena which take place at the 
same time. 

Both the chemical and the physical aspects of the various 
interfaces should be watched more carefully. Important outcomes 
of this work could be , for example, the microscopic interpre- 
tation of heat transfer coefficients at splat-chill interfaces 
on the basis of the structural features of such interfaces and 
of the capillary phenomena taking place there, or an understanding 
of the complex capillary flows which take place in shallow molten 
pools and which, in many cases, determine the final shape of the 
heat affected zone. 

Convenient starting points for the quantitative study of 
interface phenomena in RSP systems are the works by Dussan( 

1979), Timsit (1982) and Hocking (1983) , on spreading and wetting 
processes and the paper by Levich and Krylov (1969) on fluid flow 
phenomena in which surface phenomena play an important part . 

j) The success of numerical methods for the prediction of 
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cooling and solidification rates in conventional casting systems 
(e,g. Brody and Apelian(1982) and Dantzig and Berry (1984) ) , should 
be an incentive for a more widespread use of mathematical modeling 
in RS research. In particular, It would seem that considerable 
benefit can be obtained with a modest amount of effort from the 
application of the now well established method of weak solutions 
to the calculation of solidification processes during RS . 

k) The complicated fluid flow phenomena taking place in most 
RST configurations certainly warrant further study for their 
own sake. Needless to say, such phenomena are also important 
from a practical point of view since, ultimately, cooling and 
freezing rates are strongly influenced by the convective fluid 
motion taking place . It may well be that simplifications 
introduced for the sake of mathematical simplicity are such that 
important physical processes actually taking place are being 
disregarded. For this reason, more complete solutions of the fluid 
flow problem are needed. This study can be helped a great deal 
if some of the commercially available computer programs for fluid 
flow and heat transfer calculations are used. One good package 
which has produced much of the results in our group at MIT is the 
one prepared by Patankar and Spalding, In Sec(4A.l) we review 
the basic features of the Patankar -Spalding algorithm hoping that 
it finds more use within the RS community. Furthermore , in 
Sec (5. 7) we present the FORTRAN listing of a program we 
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constructed, based on the P-S algorithm, for the calculation of 
heat transfer by conduction-convection in a system where the flow 
is given by lubrication theory. It should be noted that this 
program is only a very simplified version of the P-S algorithm 
and is included here more to indicate the main ideas involved 
in the commercially available version. 

1) Last but not least is the problem of compacting the products 
of RS . Since engineering components almost always will come in 
sizes much larger than those typical of rapidly solidified samples, 
some form of compaction is unavoidable. Although a few RS 
techniques avoid this problem (e.g. liquid dynamic compaction), most 
do not. Although compacting is a very complex process, it is also 
subject to the laws of continuum mechanics. We can foresee increased 
application of both finite difference and finite element methods 
for the solution of compaction problems . 


Chapter 5 


COMPUTER PROGRAMS 


In this chapter we present the FORTRAN listings of some of 
the computer programs we have developed to perform our calculations 
of RSP systems, particularly those described in Sec(3.2)and 
(3.3). We have also included , for completeness, additional listings 
of other programs we have developed and found useful in gaining 
insight into the complex problems of RST. All the programs are 
fully commented for easier use. 

In Sec (5.1) we present a program capable of performing the 
thermal calculations for a wide range of RS configurations based 
on the assumption of Newtonian cooling conditions (Sec (3. 2)). A 
great deal of useful information can be obtained from this 
program regarding cooling and freezing rates for specific systems. 

In Secs (5. 2) and (5.3) computer programs for the calculation of 
temperature profiles and freezing rates in semiinfinite media 
under various boundary conditions, are given. The codes are simply 
programmed versions of the equations presented in Sec(3A.2). 


121 



In Sec (5.4) we present the program for the calculations of the 
PFMS system described in Sec (3. 3.1). This program is the most 
important one in the set since it produced the bulk of the results 
reported in this thesis. The program can be made to run in an 
iterative fashion by modifying a selected process parameter 
at a time, incrementally, until the overall mass balance is 
satisfied. However, we have found this process to be very 
expensive and not always convergent. Instead, we have found 
more convenient to do the iterations by performing several 
runs of the program and effecting in between reasonable changes 
in the required process parameter until our mass balance was 
satisfied. This way of running the program certainly was more 
economic. In any case it is relatively straightforward to make 
the program to perform the iterations automatically. 

Section(5.5) presents the program used to perform the 
rolling force calculations described in Sec(3.3.2). The program 
could also be used to perform calculations in related systems 
such as for the rolling of thin sheets. Section(5.6) presents 
a very simple program for the solution of systems of linear 
algebraic equations with tridiagonal matrices. This code can 
be used as the basis of an algorithm for the solution of 
transport phenomena problems in discrete form. 

Finally, in Sec (5. 7) we include a sample program we have 
constructed to perform calculations of temperatures in two- 
dimensional domains when heat is being transferred both by 
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conduction and by convection. In this program the fluid flow field 
is not calculated numerically but through closed form expressions 
obtained from lubrication theory. It is included here just to 
give an idea of the nature of one of the most widely used, 
commercially available computer programs for heat transfer and 
fluid flow calculations. 
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o o 


5.1.- Program RSPNN . Heat Transfer during RSP under 
Newtonian Cooling. 


C THIS PROGRAM COMPUTES TEMPERATURE , COOL ING RATE , SOL I F I ED FRAC- 

c TI0N A ND FREEZING RATE AS A FUNCTION OF TIME FOR VARIOUS RAPID 

c SOLIDIFICATION PROCESSING CONF I GURAT IONS . FOR SIMPLICITY, COOLING 

c Is ASSUMED TO OCCUR ACCORDING TO NEWTON ' S LAW. THE INPUT DATA 

C ARE AS FOLLOWS; TP=POURING TEMPERATURE . TO = TEMPERATURE OF THE 

C QUENCHING MEDIUM , H = HEAT TRANSFER COEFFICIENT ( CAL/CM . CM . S . C ) . 

c RH0= MELT DENSITY (G/CC) , CP =ME LT SPECIFIC HEAT (CAL/G.C) . 

c TF=MELTING POINT (C) , HF = HE AT OF FUSION (CAL/G) , 2R = SMALLEST 

C DIMENSION OF THE SPLAT EXCEPT FOR THE CASE OF THE SLAB COOLED 

C THROUGH ONE SIDE ONLY WHERE R - THICKNESS OF THE SPLAT (CM). 

c A0V = RATIO OF HEAT TRANSFER AREA TO SPECIMEN VOLUME (1/CM), 

c T I ME = TIME (S), T = TEMPERATURE (C), CR = COOLING RATE (C/S) # FS = 

c FRACTION SOLIDIFIED (-) . FR = F RE E Z I NG RATE (l/S) . RR = SOL I D I F I ED 

c DISTANCE (CM) ; GROWTH RATE (CM/S) ; GP=GEOMETR I C PARAMETER (-). 

C THE PARAMETERS FOR THE GAS ARE AS FOLLOWS: RHOG= GAS DENSITY , 

C TCG = THERMA L CONDUCTIVITY OF GAS , CPG = SPECIFIC HEAT OF GAS , 

c E MUG “VI SCOSI TY OF GAS . VG = RE LAT I VE VELOCITY OF GAS. 


TP=840. 

T0 = 0 . 

RHO = 2 . 7 
CP = 0. 206 
TF "660 . 

HF =95 . 

TCG=0. 000042 

RH0G=0. 00163 

CPG = 0 . 124 

EMUG=0. 000225 

VG= 10000.0 

COEFF = 100 . 0 

EN=0. 375 

RST AR =0 . 0025 

DR = 0 . 0005 

R I NI =R$TAR + DR 

WR I TE ( 6 , 50 ) RHO.CP.TF.HF 

WR I TE ( 6 , 50 ) TP. TO 

WR I TE { 6 , 50 ) RHOG.CPG, TCG, EMUG 

WR I TE ( 6 , 50 ) VG.DR.RINI 

WR I TE ( 6 , 50 ) COEFF. EN 

WR I TE ( 6 , 500 ) 

R=RSTAR 
DO 7 J= 1 , 21 
R= R + DR 

RE =2. *R*RHOG+ VG/EMUG 
PR -CPG + EMUG/ TCG 

H=TCG*(2. + 0-6*(RE+* .5)*(PR**.333) )/(2. ♦ R ) 

AOV = 3 . /R 
GP= 1 . /3 . 

- AOV = 3 . /R (SPHERE); 2./R (CYLINDER); 1 . /R (SLAB) 
-GP=1./3. (SPHERE) ; 1./2. (CYLINDER) ; 1. (SLAB) 

DTI ME =0.0005 
WR I TE(6 , 50 ) RE, PR 
WR I T E ( 6 , 50 ) H 
WRITE(6 , 50) R.AOV.GP 
TIME =0 . 0 
DO 1 1=1 . 1O00 


RSP00010 

RSP00020 

RSP00030 

RSP00040 

RSP00050 

RSP00060 

RSP00070 

RSP00080 

RSP00090 

RSP00100 

RSP001 10 

RSP00120 

RSP00130 

RSP00140 

RSP00150 

RSPOO160 

RSP00170 

RSP00180 

RSP00190 

RSP00200 

RSP002 10 

RSP00220 

RSP00230 

RSP00240 

RSP00250 

RSP00260 

RSP00270 

RSP00280 

RSP00290 

RSP00300 

RSP003 10 

RSP00320 

RSP00330 

RSP00340 

RSP00350 

RSP00360 

RSP00370 

RSP00380 

RSP00390 

RSP00400 

RSP004 10 

RSP00420 

RSP00430 

RSP00440 

RSP00450 

RSP00460 

RSP00470 

RSP00480 

RSP00490 

RSP00500 

RSP00510 

RSP00520 

RSP00530 

RSP00540 

RSP00550 
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c 


1 

2 


C 


3 

4 


5 

6 


7 

50 

100 

200 

300 

500 


T IME =T IME + DTI ME 

RSP00560 

TA0=(RH0*CP)/( A0V*H) 

RSP00570 

T = (TP-TO)»EXP(-TIME/TAO) + TO 

RSP00580 

CR = - (TP-TO)*EXP( -TIME /T AO) /T AO 

RSP00590 

WR I TE ( 6 , lOO) TIME.T.CR 

RSP00600 

IF(T.LE.TF) GO TO 2 

RSP00610 

CONTINUE 

RSP00620 

CONTINUE 

RSP00630 

DAS=COEFF*( (ABS(CR ) )**( - EN) ) 

RSP00640 

WR I TE ( 6 , lOO) TIME.T.CR 

RSPOO650 

WR I TE ( 6 , lOO) DAS 

RSPOO660 

TSS = T I ME 

RSP00670 

WR I TE ( 6 , 300 ) TSS 

RSP00680 

FS0=0 . 0 

RSP00690 

00 3 1=1, 1O0O 

RSP00700 

T I ME =T IME + DTI ME 

RSP00710 

C1 = ( 1 . / ( RHO*HF ) ) * AO V 

RSP00720 

C2=H* ( TF - TO ) 

RSP00730 

FS = C 1 *C2 * ( f IME-TSS ) 

RSP0O740 

FR = C 1 *C2 

RSP00750 

IF(FS.GE. 0.99999) FS=0. 99999 

RSP00760 

RR= R * ( 1 . - ( 1 . -FS)**GP) 

RSP00770 

GR= <GP*R)*( 1 ./( 1 . -FS)»*( 1 . -GP) )*FR 

RSP00780 

WRI TE ( 6 , 200 ) TIME.FS.FR.RR.GR 

RSP00790 

IF(FSO . LE . 0. 5 . AND . FS . GT . 0. 5 ) WRITE(6.200) TIME.FS 

RSP00800 

IF(FS.GE. 0.99999) GO TO 4 

RSP00810 

FSO=FS 

RSP0O820 

CONTINUE 

R5P00830 

CONTINUE 

RSP00840 

TES=TIME 

RSP00850 

WRITE ( 6 , 300 ) TES 

RSP00860 

DO 5 1=1 . 1000 

RSP00870 

TIME=TIME+DTIME 

RSP00880 

T A0= ( RHO*CP )/(AOV*H) 

R5P00890 

T=(TF-TO)*EXP( - (TIME-TES )/TAO ) + TO 

RSP00900 

CR=-(TF-T0)*EXP( -( TIME -TES )/TAO )/TA0 

RSP00910 

WR I TE ( 6 , lOO ) TIME.T.CR 

RSP00920 

I F ( T . LE . 200 . ) GO TO 6 

RSP00930 

CONTINUE 

RSP00940 

CONTINUE 

R5P00950 

WR I T E ( 6 , 1 0O ) TIME.T.CR 

RSP00960 

WRITE (6. 500) 

RSP00970 

CONTINUE 

RSP00980 

FORMAT (5X.4F 16.8) 

RSP00990 

FORMAT ( 10X, 3F20.9) 

RSPOIOOO 

FORMAT (2X.5F14.7) 

RSP010 10 

FORMAT (25X.F15.8) 

RSP01020 

FORMAT ( / ) 

RSP01030 

STOP 

RSP01040 

END 

RSP01050 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


5.2.- Program NEUMANN . Neumann’s Solution to the Stefan 
Problem. 


__ r this program computes temperatures, cooling and freezing 

RATES FOR THE CLASSICAL ONE-DIMENSIONAL STEFAN PROBLEM 

ACCORDING TO NEUMANN'S SOLUTION FOR THE SEMI - I NF INTE MEDIUM. 

SEE CARSLAW AND UAEGER (1959). CH. 11 . 

THE INPUT AND OUTPUT VARIABLES ARE AS FOLLOWS: TCL = THERMAL 

CONDUCTIVITY OF LIQUID . TCS = THERMAL CONDUCTIVITY OF SOLID , 

TDIFL * THERMAL DIFFUSIVITY OF LIQUID , TDIFS = THERMAL 

-DIFFUSIVITY OF SOLID . SPHTL = SPECIFIC HEAT OF LIQUID . SPHTS 

- = SPECIFIC HEAT OF SOLID , TINI = INITIAL TEMPERATURE OF THE 

MEL T , TWALL = WALL TEMPERATURE , TF * MELTING POINT . HEATF = 

LATENT HEAT OF FUSION , X = DISTANCE , T = TIME , XS = 

SOLIDIFIED THICKNESS , TEMP = TEMPERATURE . QUANTITIES 

IN CGS UNITS . TEMPERATURES IN CELSIUS . ENERGY IN CALORIES- 

THE DATA SHOWN CORRESPOND TO ALUMINUM SOLIDIFYING AGAINST 

A COPPER CHILL . 

TCL = .5 
TCS= . 5 
TD I FL= .925 
TD I F 5= . 925 
SPHTL= . 2 
SPHTS= . 2 
T I NI =760 . 

TWALL* 25 . 

TF = 660 . 

‘ HE ATF = 95 . 

DALAMB=0 . OO 1 
ALAMB=0 . 

DO 1 1=1 . 1000 

ALAMB=ALAMB + DALAMB 

FLHS= EXP ( -ALAMB* ALAfoB )/ERF(ALAMB) 

COE F 1 = TCL *SQRT(TDIFS)/(TCS*SQRT( TDIFL) ) 

C0EF2= (TINI - TF )/ ( TF - TWALL) 

C0EF3= EXP ( “TDI FS* ALAMB* ALAMB/ TD I FL ) 

C0EF4= 1. - ERF( ALAMB *SQRT (TDIFS/TDIFL) ) 

SLHS= - COEF 1*C0EF2*C0EF3/C0EF4 

RHS=ALAMB*HE ATF * SORT ( 3 . 14 16 )/( SPHTS + ( TF -TWALL ) ) 

DIF* FLHS + SLHS - RHS 
IF(DIF.LE. 0.001 ) GO TO 2 

1 CONTINUE 

2 CONTINUE 
ALAMB - A LAMB 

WR I TE ( 6 , 3 ) ALAMB 

3 FORMAT ( 18X , F 15 . 8 ) 

T = 0 . 

DO 5 J= 1 , 10 
DT =0.001 
T = T + DT 
X *0 . 

XS = 2. *ALAMB*SQRT (TDIFS*T ) 

WR ITE ( 6,6) XS.T 
DO 4 K = 1 . 6 
DX= . 05 

TS=((TF - TWALL )/ERF ( ALAMB )) *ERF ( X/( 2 . *SQRT( TDIFS*T )) ) + TWALL 
TL= TINI - (TINI-TF)M 1.-ERF(X/(2.*SQRT(TDIFL*T))))/< 1.- 
1 ERF(ALAMB*SQRT(TDIFS/TDIFL) ) ) 

IF(TS.GE.TF) GO TO 22 
11 CONTINUE 

TEMP=TS 
GO TO 33 
22 CONTINUE 

TEMP=TL 

33 CONTINUE 

WRITE (6,6) X , TEMP 
X=X+DX 

4 CONTINUE 

5 CONTINUE 

WR I T E ( 6 , 7 ) DT.DX.J 

6 FORMAT ( 5X . F 15 . 5 , 5X , F20. 8 ) 

7 FORMAT (/10X.2F 15.8, 5X . 13/) 

STOP 

END 


NEUOOO 1 O 

NEU00020 

NEU00030 

NEU00040 

NEU00050 

NEU00060 

NEU00070 

NEU00080 

NEU00090 

NEU00100 

NEU001 10 

NEU00120 

NEUOO 1 30 

NEU00140 

NEUOO 150 

NEUOO 160 

NEUOO 170 

NEUOO 180 

NEUOO 190 

NEU00200 

NEU002 1 0 

NEU00220 

NEU00230 

NEUOO240 

NEU00250 

NEU00260 

NEU00270 

NEU00280 

NEU00290 

NEU00300 

NEU00310 

NEU00320 

NEU00330 

NEU00340 

NEU00350 

NEU00360 

NEU00370 

NEU00380 

NEU00390 

NEU00400 

NEU004 lO 

NEU00420 

NEU00430 

NEU00440 

NEU00450 

NEU00460 

NEU00470 

NEU00480 

NEU00490 

NEU00500 

NEU00510 

NEU00520 

NEU00530 

NEU00540 

NEU00550 

NEU00560 

NEU00570 

NEU0O580 

NEU00590 

NEU00600 

NEU00610 

NEU00620 

NEUOO630 

NEU00640 

NEUOO650 

NEU00660 

NEU00670 

NEU00680 

NEU00690 

NEU00700 

NEU007 10 
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5.3 


Program SCHWARZ . Schwarz f s Solution of the Stefan 
Problem. 


C THIS PROGRAM COMPUTES TEMPERATURES. COOLING AND FREEZING 

C RATES FOR THE CLASSICAL ONE -DIMENSIONAL STEFAN PROBLEM 

C ACCORDING TO SCHWARZ'S SOLUTION FOR A SEMI - INF INTE MEDIUM. 

C -SEE CARSLAW AND JAEGER (1959) , CH. 11 . 

C THE VARIABLES HERE ARE THE SAME AS IN PROGRAM NEUMANN 

C (SEC. 5. 2) EXCEPT THAT THE THERMAL CONDUCTIVITY AND DIFFUSIVITY 

c OF THE MOLD. RESPECTIVELY TCMOLD AND TDIFM , ARE ALSO INCLUDED. 

c MOREOVER, V HERE DENOTES THE VELOCITY OF THE SOLIDIFICATION 

c INTERFACE . THE UNITS ARE THE SAME AS IN NEUMANN . 


T CMOLD= . 94 
TCL = . 29 
TCS = . 53 
TDI FM= 1 . 1 
RHOS = 2.8 
RHOL = 2 . 8 
CPS= . 23 
CPL = . 26 

TD IFL = TCL/( RHOL *CPL ) 

TDIFS = TCS/ ( RHOS + CPS ) 

SPHT L =CPL 
SPHTS=CPS 
T I N I = 7 OO . 

TWALL=25 . 

TF =660 . 

HEATF = 95 . 

DALAMB=0 . OOI 
ALAMB=0 . 

DO 1 1=1 , 1000 
ALAMB = ALAMB + DALAMB 

COE 1 = TCMOLD* SORT (TDIFS ) *EXP ( - ALAMB * ALAMB ) 
C0E2=TCS*SQRT( TDIFM) + TCMOLD * SORT ( TD I FS )* ERF ( ALAMB ) 
F LHS= C0E1/C0E2 

COEF 1 = TCL * SORT (TDIFS)/ (TCS* SORT ( TD I F L ) ) 

COE F2 = (TINI - TF ) / ( TF - TWALL) 

COEF 3= EXP( -TDIFS * ALAMB* ALAMB/TDIFL) 

COEF 4= 1. - ERF(ALAMB* SORT (TDIFS/TDIFL ) ) 

SLHS= - COEF 1*COEF2*COEF3/COEF4 

RHS = ALAMB*HEATF*SORT( 3 . 1416 )/( SPHTS*( TF- TWALL ) ) 

D I F = FLHS + SLHS - RHS 
I F ( D I F . LE .0.001 ) GO TO 2 

1 CONTINUE 

2 CONTINUE 
ALAMB= ALAMB 

WR I TE ( 6 , 3 ) DIF, ALAMB 

3 FORMAT ( 18X , 2F 15 . 8 ) 

T =0 . 

DO 5 J=1 . 10 
DT=0. 01 
T=T + DT 
X=0 . 

XS = 2 . * ALAMB* SORT ( TDI FS*T) 

V= ALAMB* SORT (TDIFS/T )*TDIFS 
WR I TE ( 6 , 7 ) XS.V.T 
DO 4 K=1 , 10 
DX=0 . 05 


SCHOOO 1 O 
SCH00020 
SCH00030 
SCH00040 
SCH00050 
SCH00060 
SCH00070 
SCH00080 
SCH00090 
SCHOOIOO 
SCH001 10 
SCH00120 
SCHOO 1 30 
SCH00140 
SCHOO 150 
SCHOO 160 
SCHOO 170 
SCHOO 180 
SCHOO 190 
SCH00200 
SCH002 10 
SCH00220 
SCH00230 
SCH00240 
SCH00250 
SCH00260 
SCH00270 
SCH00280 
SCH00290 
SCH00300 
SCH00310 
SCH00320 
SCH00330 
SCH00340 
SCH00350 
SCH00360 
SCH00370 
SCH00380 
SCH00390 
SCH00400 
SCH004 lO 
SCH00420 
SCH00430 
SCH00440 
SCH00450 
SCH00460 
SCH00470 
SCH00480 
SCH00490 
SCH00500 
SCH00510 
SCH00520 
SCH00530 
SCH00540 
SCH00550 
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00 -4 CT) Ul 


1 1 


22 


33 

44 


4 


AM 1 = TCS*SQRT(TDIFM) + (TF-TWALL) 

AM2= TCS*SQRT(TDIFM) + TCMOLD* SORT (TDIFS)*ERF( A LAMB ) 

AMOLD= AM1/AM2 

TMOLD= AMOLD* ( 1 .+ERF(X/(2.+ SORT (TDIFM+T)))) + TWALL 
AS 1 = ( TF -TWALL )/ AM2 
AS2= TCS*SQRT(TDIFM) 

AS3= T CMOLD *SQRT ( TD I FS ) 

TS= TWALL + AS 1 * ( AS2 + AS3+ERF ( X/ ( 2 . *SORT( TDI FS*T ) ) ) ) 

AL1 = (TINI-TF)/( 1 . -ERF(ALAMB* SORT (TDIFS/TOIFL) ) ) 

TL= TINI - AL 1 *( 1 . -ERF(X/(2 . * SORT ( TD I F L * T ) ) ) ) 

IF(X . LE .0 . ) GO TO 33 

I F (TS . GE . TF ) GO TO 22 

CONTINUE 

TEMP = T S 

GO TO 44 

CONTINUE 

TEMP = T L 

GO TO 44 

CONTINUE 

TEMP=TMOLD 

CONTINUE 

WRITE (6,6) X , TEMP 

X=X+DX 

CONTINUE 

AL2=EXP( - ALAMB + ALAMB * TD I FS/TDIFL ) 

AL3 = 4 . * ALAMB + ALAMB * ALAMB * TD I F S *SORT (TDIFS/TDIFL )/ SORT (3. 14 1592 ) 
AL=-AL 1*AL2*AL3 
CR=AL/ ( XS*XS ) 

WRITE(6,8) CR 
CONTINUE 

FORMAT ( 5X , F 1 5 . 5 , 5X . F 20 . 8 ) 

FORMAT ( 1X.3F20-9) 

FORMAT ( 13X, F25.9) 

STOP 

END 


SCH00560 
SCH00570 
SCH00580 
SCH00590 
SCH00600 
SCH0O610 
SCH00620 
SCHO0630 
SCH00640 
SCH00G50 
SCH00660 
SCH00670 
SCH00680 
SCH00690 
SCH0O700 
SCH007 10 
SCH00720 
SCH00730 
SCH00740 
SCH00750 
SCH00760 
SCH00770 
SCH00780 
SCH00790 
SCH00800 
SCH008 10 
SCH00820 
SCH00830 
SCH00840 
SCH00850 
SCH00860 
SCH00870 
SCH00880 
SCH00890 
SCH00900 
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5.4.- Program PFMS . Calculation of Heat Transfer and Fluid 


Flow in the Planar Flow Melt Spinning Process. 


C- 

c- 

C- 

c- 

c- 

c- 

C- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 


c- 

c- 

c- 

c- 

c- 

c- 


c 

c 

c 


-THIS PROGRAM COMPUTES SOLIDIFIED THICKNESS FOR THE CASE OF SINGLE CCS00010 
-ROLL STRIP CASTING OF METALS. THE ALGORITHM IS BASED ON THE ENTHALP YCCS00020 
-METHOD AND THE PARTICULAR SCHEME USED HAS BEEN AN EXPLICIT ONE. CCS00030 

-ENTHALPIES (AND THUS TEMPERATURES) ARE COMPUTED EXPLICITLY FOR A CCS00040 

-ROW OF GRID POINTS AT A GIVEN DOWNSTREAM LOCATION BY USING THE VALUE CCS00050 
-OF ENTHALPY AND TEMPERATURE OF THE GRID POINTS ALONG THE INMEDI ATLY CCS00060 
-PRECEEDING UPSTREAM LOCATION. CCS00070 

-HERE WE ASSUME "PLUG FLOW" TYPE MOTION OF THE METAL BEING CAST. CCS00080 

-HOWEVER, VELOCITY PROFILES ARE COMPUTED AT EVERY DOWNSTREAM LO- CCS00090 

-CATION TAKING INTO ACCOUNT THE PRESENCE OF A SOLIDIFIED SHELL. CCSOOIOO 

-DURING THE LAST STAGES OF SOLIDIFICATION, FLUID FLOW AND HEAT CCSOOIIO 

-TRANSFER INTERACT MUCH STRONGER SINCE THE FREE MELT-GAS SURFACE CCS00120 

-IS COMPUTED FROM THE SOLIDIFIED THICKNESS AND THIS, IN TURN . CCS00130 

-DEPENDS ON THE PRECISE LOCATION OF THE FREE BOUNDARY. CCS00140 

REAL UO( 32) ,UN(32 ) ,HO( 32) ,HN( 32) ,CTU(32 ) .TEMP ( 32) , TEMPO ( 32 ) CCS00150 

REAL TINI (32 ) ,CR( 32 ) . VX( 32 ) . VXO(32 ) , STRE AM( 32 ) CCS00160 

--THE VARIABLES IN THE ARRAYS ARE: UO =OLD PSEUDO TEMPERATURE (SEE CCS00170 
--BELOW); UN = NEW PSEUDO TEMPERATURE; HO = OLD ENTHALPY ; HN = NEW CCS00180 

--ENTHALPY; CTU = C TIMES U (SEE BELOW); TEMP = ACTUAL TEMPERATURE; CCS00190 

--TEMPO = OLD TEMPERATURE; TINI = INITIAL TEMPERATURE ; CR = COOLING CCS00200 
--RATE ; VX * X-COMPONENT OF VELOCITY ; VXO = OLD VELOCITY ; STREAM CCS00210 

-- = STREAM FUNCTION ( D I MENT I ONLE SS ) . CCS00220 

TFUS =1325.0 CCS00230 

HFUS = 7 1 . 7 CCS00240 

TC=0 . 07 1 7 CCS00250 

CP = 0. 1434 CCS00260 

RH0=8 . 5 CCS00270 

EMU=0 . 046 CCS00280 

SIGMA =1778.0 . CCS00290 

TCR =0 . 16 CCS00300 

TCR = 0 . 93 CCS00310 

RH0R=7 . 86 CCS00320 

RHOR x 8 . 94 CCS00330 

CPR =0 . 1 5 CCS00340 

CPR=0 . 09 1 4 CCS00350 

TDR= TCR/ ( RHOR* CPR ) CCS00360 

WR ITE(6, 1 10) TFUS .HFUS.TC , CP CCS00370 

WRI TE ( 6 , 110) RHO, EMU. SIGMA CCS0O380 

WR I TE ( 6 . 1 10) TCR, RHOR. CPR. TOR CCS00396 

-- QUANTITIES WILL BE GIVEN IN CGS UNITS. CCS00400 

--ENERGIES WILL BE GIVEN IN CALORIES AND TEMPERATURES IN CELSIUS. CCS00410 

--THE PHYSICAL PROPERTY DATA ARE AS FOLLOWS: TFUS = MELTING POINT OF CCS00420 
--THE SUBSTANCE BEING CAST (NOTE; IN THE CASE OF AN ALLOY, INSTEAD OF CCS00430 
--TFUS A LIQUIDUS AND A SOLIDUS HAVE TO BE SPECIFIED); HFUS = LATENT CCS00440 
--HEAT OF FUSION ; TC = THERMAL CONDUCTIVITY ; CP = SPECIFIC HEAT ; CCS00450 

-- RHO = DENSITY ; EMU = VISCOSITY ; SIGMA = SURFACE TENSION ; TCR = CCS00460 

-- THERMAL CONDUCTIVITY OF ROLL ; RHOR = DENSITY OF ROLL ; CPR = SPE- CCS00470 
--CIFIC HEAT OF ROLL ; TDR = THERMAL DIFFUSIVITY OF ROLL . CCS00480 

PI *3 . 14159 CCS00490 

TWALL = 25 . 0 CCS00500 

TBULK= 1440.0 CCS00510 

BN=0 . 064 CCS00520 

HI N=0 . 0300 CCS00530 

W=0 . 635 CCS00540 

OMEGA = 1200.0 CCS00550 
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RR= 1 2 . 7 
0 = 3.86 
PL*0. 29 
THETA*20 . 0 

VRXR= 0MEGA*2 . 0*P I *RR/60 . 0 
QPUW= Q/W 

COEF I = 3 . 0*EMU/SIGMA 
TW = TWALL*TCR 
TF = TFUS*TC 
TB = TBULK*TC 
HF= HFUS*RH0 
XI=-0. 60 

BN 
PL 


XB= XI + 
X0 = XI + 
XF =0 . 0 
X INI = XI 
XBRE = XB 
XDET= XO 
XFIN = XF 


XI 

XI 

XI 

XI 


F0= HF«Q/(W*( 1 .5*PL) ) 

TGR = FO/TCR 
H 10= HIN 

H20= - T AN( THETA*2.0*PI/ 360 . 0 ) 
H30= 1 .0/ ( PL/2 . ) 

SLI PE =0 . 25 
HTCO= 1 .035 
YFS = 0 . 0 
P0=O.O 

PF = 10 1300000.0 
ALF A =0 . 0 
BETA=0 . 333 
GAMA = 1 .002 
WRI TE ( 6 , 1 10) 

WRITE ( 6 , 1 10) BN. HIN. W 
WRITE ( 6 . 110) OMEGA, RR,0 
WR ITE(6. 1 10) PL, THETA 
WRITE(6.110) VRXR , OPUW , COEF I 
WRITE(6, 1 10) XI.XB.XD.XF 

XINI ,XBRE .XDET .XFIN 

FO . TGR 
H 1 0 , H20 , H30 
SLIPE .HTCO 
BETA .GAMA 


TWALL . TBULK 


WRI TE ( 6 , 1 10) 

WRITE (6 . 1 10) 

WRITE (6 . 1 10) 

WRITE (6 . 1 10) 

WR I TE ( 6 , 1 10) 

C---THE PROCESS PARAMETERS ARE : TWALL = TEMPERATURE OF THE ROLL (SEE 

C BELOW ) ' TBULK = POURING TEMPERATURE ; BN = NOZZLE BREADTH : HIN = 

c WHEEL-NOZZLE GAP ; W = NOZZLE (STRIP) WIOTH ; OMEGA = ANGULAR VELO 

r CITY OF ROLL (RPM); RR = ROLL RADIUS : O = MELT FLOW RATE ; PL = 

r PUDDLE LENGTH • THETA = CONTACT ANGLE MELT-NOZZLE (DOWNSTREAM) : 

S — VRXR - SURFACE VELOCITY OF ROLL ; XINI . XBRE . XDET AND XFIN ARE . RESPECCCS0 1040 
C — TIVELY THE LOCATIONS OF : THE UPSTREAM EDGE OF PUDDLE . THE END OF CCS01050 

r NOZZLE BREADTH , THE DETACHMENT POINT AND THE (APPROXIMATE) END' 

C PO I NT OF SOLIDIFICATION ; H10.H20 AND H30 ARE USED TO INITIATE THE 

C NUMERICAL SOLUTION FOR THE FREE SURFACE ; SLIPE = SLIP EXPONENT ( 

q SEE BELOW - -LBL . 16) ; HTCO = HEAT TRANSFER COEFFICIENT J^SEE BELOW 

C , -LBL . 33) 


CCS00560 
CCS00570 
CCS00580 
CCS00590 
CCS00600 
CCS00610 
CCS00620 
CCS00630 
CCS00640 
CCS00650 
CCS00660 
CCS00670 
CCS00680 
CCS00690 
CCS00700 
CCS007 10 
CCS00720 
CCS00730 
CCS00740 
CCS00750 
CCS00760 
CCS00770 
CCS00780 
CCS00790 
CCS00800 
CCS008 10 
CCS00820 
CCS00830 
CCS00840 
CCS00850 
CCS00860 
CCS00870 
CCS00880 
CCS0O89O 
CCS00900 
CCS00910 
CCS00920 
CCS00930 
CCS00940 
CCS00950 
CCS00960 
CCS00970 
CCS00980 
CCS00990 
CC SO 1000 
- CCS01010 
CCS01020 
CCS01030 


CCS0106O 

CCS01070 

CCS01080 

CCS01090 


BETA = FRACTION TO DEFINE THE LOCATION OF SOLID-LIQUID CCSOIIOO 
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C INTERFACE (SEE AFTER LBL . 37) ; GAMA = COEFFICIENT TO DECIDE IF A 

C GIVEN GRID POINT HAS SOLIDIFIED (SEE AFTER LBL. 37 BELOW). 

C 

N= 31 
NMAX=N 
NM I N 1 = N- 1 
NT * 8000 
C 

c THE GRID PARAMETERS ARE: N = NUMBER OF POINTS ALONG THE Y (VER- 

C TICAL, PERPENDICULAR TO THE WHEEL) DIRECTION , NT = NUMBER OF 

C POINTS ALONG THE X (DOWNSTREAM) DIRECTION . 

C THE FOLLOWING LOOP INITIALIZES BOTH TEMPERATURE AND ENTHALPY. 

DO 10 1= 1 . N 
U0( I )= TB 

HO ( I ) - ( ( RHO*CP )/TC)*(UO(I)-TF) + HF 
UN( I ) s TB 
T INI ( I )= U0( I )/TC 
10 CONTINUE 

WR I TE ( 6 , 1 00 ) (TINI(I).I* 1,N) 

T I ME = O. 

X = XI 

D I STX= 0.0 

KOUNT 1=0 
KOUNT 2=0 

C THE FOLLOWING LOOP ADVANCES THE CALCULATION ALONG THE DOWNSTREAM 

C DIRECTION, JUMPING FROM A LINE OF GRID POINTS AT CONSTANT X TO THE 

C NEXT . THIS IS THE MAIN "OUTER LOOP" OF THE CALCULATION . 

C NOTE; TIME, KOUNT 1 AND KOUNT 2 ARE ONLY COUNTERS. 

DO 50 K= 1 , NT 

KOUNT 1=K0UNT 1+1 
KOUNT2=KOUNT2+ 1 

c "OLD" TEMPERATURE IS CALCULATED. 

DO 12 J= 1 , NMAX 
TEMPO ( J ) -UO ( J ) /TC 
12 CONTINUE 

C IN THE FOLLOWING THE GRID SPACING IS COMPUTED, SINCE THE METHOD IS 

C AN EXPLICIT ONE , DX AND DY ARE RELATED BY THE STABILITY CONDI- 

C TION. THIS PART OF THE PROGRAM (UP TO LBL. 15) COMPUTES THE LOCA- 

C TION OF THE UPPER BOUNDARY. THEN . UP TO LBL 16. THE GEOMETRICAL 

C PARAMETERS OF THE GRID ARE SET. 

C THE GAS-MELT INTERFACE HAS TO BE SPECIFIED ( DI STX . GT . XDET ) EITHER 

C BY SOME (ASSUMED) SHAPE OR BY NUMERICALLY SOLVING THE CAPILLARY 

C EQUATIONS (PREFERRED). 

C THE NEXT THREE CONDITIONS CONTROL WHERE THE MELT-GAS INTERFACE HAS 

C TO BE CALCULATED DEPENDING ON THE PROGRESS OF SOLIDIFICATION AND 

C ON IF D I STX IS .GT. XDET . 

H = HIN 

HMYFS = H - YFS 
IF(HMYFS.LE. 0.0001) GO TO 51 
I F ( UN( NMAX ) . LE . GAMA *TF ) H = HF IN 
I F ( UN( NMAX ) . LE . GAMA * TF ) GO TO 15 
IF(DISTX. LT .XDET) GO TO 15 
14 CONTINUE 

HMYFSC = HMYF S*HMYFS + HMYF S 
HMY F 5S = HMY F S * HMYFS 


CCS01 1 10 
CCS01 120 
CCS01 130 
CCS01 140 
CCS01 150 
CCS01 160 
CCS01 170 
CCS01 180 
CCS01 190 
CCS01200 
CCS01210 
CCS01220 
CCS01230 
CCS01240 
CCS01250 
CCS01260 
CCS01270 
CCS01280 
CCS01290 
CCS01300 
CCS01310 
CCS01320 
CCS01330 
CCS01340 
CCS01350 
CCS01360 
CCS01370 
CCS01380 
CCS01390 
CCS01400 
CCS01410 
CCS01420 
CCS01430 
CCS01440 
CCS01450 
CCS01460 
CCS01470 
CCS01480 
CCS01490 
CCS01500 
CCS01510 
CCS01520 
CCS01530 
CCS01540 
CCS01550 
CCS01560 
CCS01570 
CCS01580 
CCS01590 
CCS01600 
CCS01610 
CCS01620 
CCS01630 
CCS01640 
CCS01650 
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FUNCT = COEF I * ( ( QLPUW/HMY FSC ) - ( VRX/HMY FSS ) ) 


CCS01660 


H3 = H30 + FUNCT *DX 


CCSO 1670 


H2= H20 + H3*DX 


CCS01680 


H 1 = H10 + H2*DX 


CCSO 1690 


IF (H2.GE .0.0001 ) GO TO 51 


CCSO 1 700 


H10=H1 


CCS017 10 


H20=H2 


CCS0 1 7 20 


H30=H3 


CCSO 1 7 30 


HMYFS = HI 


CCS0 1 7 40 


H= HMYFS + YFS 


CCSO 1 750 


IF(H.GT.HIN) H = HIN 


CCSO 1 760 

15 

CONTINUE 


CCSO 1 770 


VELAV= Q/(W*HIN) 


CCSO 1 7 80 


V=Q/(W*H) 


CCSO 1 790 


DY=H IN/FLOAT ( NMIN1 ) 


CCSO 1 800 


DX= (DY *DY ) ‘VELAV/2 . 


CCS018 10 


DT=DX/VELAV 


CCSO 1820 


C-DT / ( DY *DY ) 


CCSO 1830 


T I ME = T I ME +DT 


CCSO 1840 


DISTX=DISTX+DX 


CCSO 1850 


X=X + OX 


CCSO 1860 

16 

CONTINUE 


CCSO 1870 

c 



CCSO 1 880 


SL 1 PC= OISTX/XFIN 


CCSO 1890 


IF(DISTX.GE XFIN) SLIPC 

= 1.00 

CCSO 1 900 


VRX* VRXR*SLIPC**SLIPE 


CCS01910 

c 



CCSO 1920 

C THE 

NEXT 4 STATEMENTS SERVE TO ACCOUNT FOR THE PRESENCE 

OF 

CCSO 1930 

C THE 

FREE SURFACE AFTER DETACHMENT. GRID POINTS OUTSIDE 

THE 

CCSO 1940 

C- - -MELT 

ARE NOT COMPUTED THERE. 


CCSO 1950 


DI F = ( HIN + 0 . 5*DY - H )/DY 


CCSO 1 960 


KDIF=INT(DIF ) 


CCSO 1970 


NMAXM 1 = ( N- 1 ) - KDIF 


CCSO 1 980 


NMAX= NMAXM 1 + 1 


CCSO 1 990 


NMAXM2 = NMAXM 1 - 1 


CCS02000 


Y =0 . 0 


CCS020 1 0 

C WITH 

1 THIS LOOP WE ADVANCE THE CALCULATION IN THE VERTICAL (Y) 

CCS02020 

C DIRECTION. FROM THE EON. V+DH/DX = K D(DT/DY)/DY WE CONSTRUCT 

CCS02030 

C THE 

EXPLICIT F.D. EON.; HN = HO + C + ( U0( J - 1 ) * 2U0 ( J ) +U0 ( J + 1 ) ) . 

CCS02040 

C THE 

RESULT OF THIS CALCULATION IS HN(K) . 


CCS02050 

C THIS 

( is THE MAIN ‘'INNER" LOOP OF THE PROGRAM. 


CCS02060 


DO 20 1=2, NMAXM 1 


CCS02070 


IP1 = I + 1 


CCS02080 


IM 1 3 I - 1 


CCS02090 


CTU ( I )=C + (UO( IM1 )+U0( IP1 )-2 . *U0( I ) ) 


CCS02 lOO 


HN( I ) =H0 ( I ) +C TU ( I ) 


CCS02 1 lO 

20 

CONTINUE 


CCS02 120 

C THE 

FOLLOWING STATEMENTS (UP TO 30) SIMPLY COMPUTE THE 

PSEUDOTEM- 

CCS02130 

C PERATURES CORRESPONDING TO THE JUST CALCULATED VALUES OF HN BY 

CCS02 140 

C USING THE THERMODYNAMIC RELATIONSHIP BETWEEN ENTHALPY 

AND TEMPE- 

CCS02150 

C RATURE . 


CCS02160 


00 30 I =2 , NMAXM1 


CCS02170 


IF (HN( I ) . GT . HF )G0 TO 27 


CCS02180 

23 

CONTINUE 


CCS02 190 


IF ( HN( I).LE.HF.AND. HN( I ) . GE . 0 . )G0 TO 26 


CCS02200 
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o o 


w 

24 

CONTINUE 

CCS022 10 



UN(I)=TF + HN(I)*(TC/(RHO*CP)) 

CCS022 20 



GO TO 30 

CCS02230 


26 

CONTINUE 

CCS02240 



UN( I ) =TF 

CCS02250 

* 


GO TO 30 

CCS02260 


27 

CONTINUE 

CCS02270 



UN ( I ) = T F + ( HN( I ) -HF ) + (TC/( RHO *CP ) ) 

CCS02280 


30 

CONTINUE 

CCS02290 


C TO 

COMPUTE THE TEMPERATURE AT THE WHEEL-SPLAT INTERFACE IMPERFECT 

CCS02300 


C THERMAL CONTACT IS ASSUMED . THE TEMPERATURE AT THE SURFACE OF THE 

C WHEEL CAN BE GIVEN BY THE FORMULA FOR PEAK TEMPERATURE IN A 

C THICK SOLID UNDER A MOVING HEAT SOURCE. TWO ALTERNATIVE B.C'S FOR 

C THIS BOUNDARY ARE THAT UW = TW AND THAT UN( 1 ) = UW (VARIABLE) 

C , WHICH WOULD CORRESPOND TO IDEAL COOLING (I.E. HTCO VERY LARGE). 

C MOREOVER , HTC CAN BE MADE VARIABLE AS A FUNCTION OF VRX . 
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CCS02310 
CCS02320 
CCS02330 
CCS02340 
CCS02350 
CCS02360 
CCS02370 
CCS02380 
CCS02390 
CCS02400 
CCS024 10 
CCS02420 
CCS024 30 
CCS02440 
CCS02450 
CCS02460 
CCS02470 
CCS02480 
CCS02490 
CCS02500 
CCS025 1 0 

C CCS02520 

C THE FOLLOWING STATEMENTS (UP TO 37) INTRODUCE THE BOUNDARY COND I T I 0NCCS02530 

C ON THE NOZZLE (FREE SURFACE) SIDE OF THE GRID. NOTE THAT BENEATH CCS02540 

C THE NOZZLE BREADTH THE TEMPERATURE IS TAKEN AS THE POURING VALUE CCS02550 

C WHI LE AFTER THE BREADTH AND ALONG THE FREE SURFACE A ZERO HEAT FLOW CCS02560 


CONTINUE 

DU= ( 2 . * FO ) +SQRT ( TDR * ( D I STX/ VRXR ) /P I ) 

UWR= TW + DU 
TSR= UWR/TCR 
UW-UWR 
UW= TW 

TSURF = UW/TCR 
CHTC= 1 . 0 

I F ( UN( 1 ) . GT . GAMA *TF ) CHTC = VRX/ VRXR 
HTC = HTCO *CHTC 

UN(1) = ( (TC/DY)*UN(2) + HTC *UW ) / ( ( TC/DY ) + HTC) 
UN( 1 ) = UW+TC/TCR 

HN( 1 ) = ( (RH0 + CP)/TC)*(UN(-1) - TF) 

IF(UN( 1 ) .GE.TF) HN( 1 ) * ( ( RHO*CP )/TC ) * ( UN( 1 ) - TF ) 


HF 


C COND I T I ON IS USED. 

UN( NMAX ) =UN( NMAXM 1 ) 

HN ( NMAX ) -HN( NMAXM 1 ) 

IF(DISTX. LE . XBRE ) GO TO 36 
GO TO 37 
CONTINUE 
UN( NMAX ) = TB 

HN ( NMAX ) = ( (RHO+CP )/TC ) ♦ ( TB - TF ) 
CONTINUE 


35 

36 


HF 


37 


CCS02570 
CCS02580 
CCS02590 
CCS02600 
CCS026 10 
CCS02620 
CCS02630 
CCS02640 
CCS02650 


C THE FOLLOWING PORTION (FROM 37 UP TO 38) COMPUTES VELOCITY PROFILES CCS02660 

C IN THE SPLAT UNDER THE ASSUMPTIONS OF LUBRICATION THEORY. ACCOUNT CCS02670 

C---IS TAKEN OF THE THICKNESS SOLIDIFIED BY MEANS OF ALFA . CCS02680 

YFS=0 . 0 CCS02690 

I F ( UN( 1 ) . LE . GAMA *TF ) YFS = BETA +DY CCS02700 

DO 301 J = 2 , NMAXM 1 CCS02710 

JM 1 = J - 1 CCS02720 

JP1=J+1 CCS02730 

IF(UN(J).LE.GAMA*TF) YFS= ( FLOAT( JM 1 ) -BETA ) *DY CCS02740 
I F ( UN( J ) . LE . TF / GAMA ) YFS= ( F LOAT ( JM 1 ) +0 . 5 *BE T A ) *DYCCS02750 
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UNPR = (UN(UPt) + UN(J))/2. 

I F ( UNPR . LE . GAMA * TF ) YFS= (FLOAT ( JM 1 )+BET A ) *DY 

I F ( UNPR . LE . TF/GAMA ) YFS= ( FLOAT ( JM 1 )+ 1 . 5*BET A ) *DY 

IF ( UN( J ) . GT . GAMA*TF ) GO TO 302 

CONTINUE 

CONTINUE 

I F ( UN( NMAX ) .LE . GAMA *TF ) Y F S = ( F LOAT ( NMAXM t ) - BETA ) *DY 
ALFA= YFS/H 

I F ( ALF A . GE . 1.000) ALFA =0.9999 

FOA= 1. - 3 . * ALFA + 3. * ALF A * ALF A - ALF A * ALF A * ALF A 

IF(OISTX.GE.XDET) GO TO 305 
CONTINUE 

PP02E= 6.*( VRX * ( H/ 2 . ) * ( 1 . “ ALF A ) + VRXR*H*ALFA - QPUW )/ 

( ( H*H* H ) * FOA ) 

A2= PP02E 

A1= VRX/(H*( ALFA- 1 . ) ) - PP02E*H*(ALFA+1. ) 

A0= PP02E *H*H*ALFA - VRX/ ( ALFA - 1 . ) 

QLPUW = - PP02E*((H*H*H)/6.)*F0A + VRX+(H/2.)*(1.-ALFA) 

GO TO 307 

CONTINUE , . 

PP02E = ( 3 . /2 . ) ♦ ( VRX + H+ ( 1 . -ALFA ) + VRXR + H * ALF A - QPUW )/ 

( ( H*H*H ) + FOA ) 

A2= PP02E 
A 1 = - PP02E * 2 . *H 

A0= VRX - PP02E *H*H* ( ALFA + ALFA - 2 . * ALF A ) 

OLPUW= - PP02E*(2./3. )*(H*H*H)*FOA + VRX*H* ( 1 . - ALFA ) 

CONTINUE 


♦ALFA) 


PP02E*(2./3. )*(H*H*H)*FOA + VRX*H* ( 1 . - ALFA ) 


ATH= ALFA*H 
ATHS= ALFA*H* ALFA*H 
ATHC- ALFA*H*ALFA*H*ALFA*H 
PP= 2 . ♦ EMU * PP02E 
PN= PO + PP ♦DX 

IF ( PN . GE . PF . AND .PO.LT.PF) XPRESS=DISTX 
PO = PN 

Y = 0 . 0 

DO 38 1*1 .NMAX 

IF ( ALFA*H. GT .0.0. AND . Y . LT . ALFA*H) GO TO 330 

310 CONTINUE 

VX(I)= A2 * ( Y * Y ) + A 1 + Y + AO 

STREAM( I )=(VRXR*ALFA*H + (A2/3.)*( Y*Y + Y - ATHC ) + 

1 (A1/2.)+< Y*Y * ATHS ) + AO* ( Y - ATH ))/QPUW 

GO TO 350 

330 CONTINUE 

VX ( I ) = VRXR 

STREAM ( I ) = VRXR + Y / QPUW 
350 CONTINUE 

Y=Y+DY 

38 CONTINUE 

IF(UN( 1) .LE.TF.AND.UOt 1 ) .GT.TF) XSS=DISTX 
I F ( VX ( NMAX ) . GE .0.00. AND . VXO(NMAX) . LT .0.00) XSTAG2-DISTX 

I F(UN( NMAX ). LE . GAMA *TF . AND . U0( NMAX ). GT . GAMA*TF ) HFIN = H 
I F (UN( NMAX ). LE . TF/GAMA . AND . U0( NMAX ). GT . TF/GAMA ) GO TO 51 

q JO AVOID HAVING TO STORE THE WHOLE GRID, THE FOLLOWING LOOP RESETS 

C---UO TO BE THE FRESHLY CALCULATED VALUE I.E. UN . SAME WITH HN . 

39 CONTINUE 
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DO 40 1*1, NMAX CCS03310 

U0( I ) =UN( I ) CCS03320 

HO ( I )*HN( I ) CCS03330 

VXO( I ) = VX ( I ) CCS03340 

40 CONTINUE CCS03350 

C THE FOLLOWING LOOP RECOVERS THE ACTUAL TEMPERATURES AND ALSO CCS03360 

C COMPUTES THE VALUES OF THE COOLING RATE FOR EVERY GRID POINT. CCS03370 

C BESIDES, AVERAGE COOLING RATES ACROSS THE SPLAT FOR FIXED X CCS03380 

C LOCATION ARE CALCULATED. CCS03390 

DO 41 1=1 .NMAX CCS03400 

, TEMP( I )=UN( I )/TC CCS03410 

41 CONTINUE CCS03420 

SUM=0 . O CCS03430 

DO 42 1=1, NMAX CCS03440 

CR ( I )=(TEMP( I )- TEMPO ( I ) ) * ( VX ( I )/DX ) CCS03450 

SUM= SUM + CR( I ) CCS03460 

42 CONTINUE CCS03470 

AVERCR= SUM/NMAX CCS03480 

C WITH THE FOLLOWING STATEMENTS IT IS POSSIBLE TO PRODUCE OUTPUT WITH CCS03490 

C DIFFERENT FREQUENCIES IN THE DIFFERENT REGIONS OF THE DOMAIN. CCS03500 



IF(ALFA.GT.O. 20 ) GO TO 48 

CCS035 10 

44 

IF ( KOUNT 1 . NE . 200) GO TO 50 

CCS03520 

46 

GO TO 49 

CCS03530 

48 

I F ( K0UNT2 . NE . 200 ) GO TO 50 

CCS03540 

49 

CONTINUE 

CCS03550 

C444 

GO TO 50 

CCS03560 

C---I 

FINALLY RESULTS ARE WRITEN 

CCS03570 


WRITE (6, 102) K , I P 1 , IM 1 , NMAX , NMAXM 1 

CCS03580 


WRI TE ( 6 , 99 ) H.DY.DISTX.X, DX .TIME 

CCS03590 

C 

WR I TE ( 6 , lOO ) (HN(I ), 1=1 .NMAX) 

CCS03600 


WRITE(6, 123) TSURF .TSR.YFS.ALFA 

CCS036 1 0 


WRI TE ( 6 , 126) AVERCR 

CCS03620 


WRI TE ( 6 , 126) PP , PN 

CCS03630 


WR I T E ( 6 , 1 00 ) (TEMPI I ), 1= 1 .NMAX) 

CCS03640 


WR I TE ( 6 , 105 ) (CR( I ) . 1 = 1 .NMAX) 

CCS03650 


WR I TE ( 6 . lOO ) (VX( I ) . 1 = 1 .NMAX ) 

CCS03660 


WR I TE ( 6 , 107) ( STREAM* I). 1 = 1 .NMAX) 

CCS03670 


KOUNT 1=0 

CCS03680 


KOUNT 2 = 0 

CCS03690 

50 

CONTINUE 

CCS03700 

51 

CONTINUE 

CCS037 10 


OF = VRXR*W*H 

CCS03720 


ERRORQ= (0 - OF )/Q 

CCS03730 


WRITE (6. 102) K . I P 1 , I M 1 . NMAX . NMAXM 1 

CCS03740 


WR I TE ( 6 , 99 ) H.DY.DISTX.X.DX.TIME 

CCS03750 


WR I TE ( 6 , 124) TSURF. XSS . XSTAG2 . XDET 

CCS03760 


WR I TE ( 6 , 125 ) YFS.ALFA.OF.ERRORO 

CCS03770 


WRITE(6, 126) PP.PN.XPRESS. HF I N 

CCS03780 


WRI TE ( 6 , 100 ) (TEMP(I), 1=1, NMAX) 

CCS03790 


WRITE(6. 105) (CR(I). 1=1, NMAX) 

CCS03800 


WR ITE ( 6 , 100) ( VX ( I ) . I = 1 , NMAX ) 

CCS03810 


WRI TE ( 6 . 107 ) (STREAM* I ). 1 = 1 .NMAX) 

CCS03820 

99 

FORMAT ( 5X , 6F 10.6) 

CCS03830 

100 

FORMAT ( IX ,4F 15. 5) 

CCS03840 

102 

FORMAT ( 15X.5I 10) 

CCS03850 


105 FORMAT ( 1X.4E15.5) 

107 FORMAT ( 1X.4F15.8) 

1 10 FORMAT ( 5X , 4F 16.6) 

123 FORMAT ( 10X.4F15.5) 

124 FORMAT ( 5X , 5F 13.5) 

125 FORMAT ( 12X . 4F 15 . 5 ) 

126 FORMAT ( 12X , 4E 15 . 5 ) 

101 CONTINUE 

STOP 

END 


CCS03860 

CCS03870 

CCS03880 

CCS03890 

CCS03900 

CCS03910 

CCS03920 

CCS03930 

CCS03940 

CCS03950 
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5 . 5 .- 


Program PRESSTRQ . Calculation of Rolling Forces in the 
Gap of a Twin Roll RS Device . 


THIS PROGRAM COMPUTES THE PRESSURE DISTRIBUTION IN THE 

-TWIN-ROLL-GUENCHING MACHINE FOR SPLAT COOLING. THE PROGRAM 
-USES A SINGLE-STEP , FOURTH-ORDER RUNGE KUTTA METHOD AS 
-DESCRIBED IN FERZIGER J. NUMERICAL METHODS FOR ENGINEERING, 
-WILEY 1981 P 

THE VARIABLES NEEDED FOR THE CALCULATION ARE; 0. THE 

-VOLUME FLOW RATE PER UNIT WIDTH, R, THE ROLL RADIUS;HO. 

-THE MINIMUM GAP BETWEEN THE ROLLS ; OMEG , THE ANGULAR VELOCI- 
TY OF THE ROLLS; XI THE LIFT-OFF P0INT;X4,THE POINT OF INITIAL 
-CONTACT WITH THE ROLLS;XN, THE SLIP COEFF ICI ENT ; PL AND B 
-ARE THE RHEOLOGICAL PROPERTIES OF THE SPLAT. 

---PNEW IS THE REQUIRED PRESSURE AND DDX IS THE STEP SIZE USED 
-IN THE CALCULATIONS . DPDX IS THE PRESSURE GRADIENT AND X IS THE 
-SPACE AXIS IN THE ROLLING DIRECTION. 

NOTE THAT THE CALCULATED DPDX VALUES ARE BASED ON THE 

-ASSUMPTIONS OF THE THEORY OF LUBRICATION. 


DIMENSION H(3) , SL( 3 ) ,VRX(3) , VRXP( 3),A(3),C1(3) ,DPDX( 3 ) 

Q= . 97 
R=5.0 
HO =0.005 
OMEG= 160.00 

C CALCULATION OF THE LIFT-OFF POINT 

COEO= 0*60 . /( 3 . 14 16*0MEG) 

COE 1 = HO + R 

C0E2* SORT ( COE 1**2. - COEO ) 

C0E3= (COE 1 + C0E2)/2. 

X1= SORT ( ABS ( R* *2 . - C0E3**2.)) 

X4= -0 . 9 

H4= HO+R-SQRT (R*R-X4*X4) 

VIN = 0/(2. *H4 ) 

XN=0 . O 
PL = 5 . 0 

C THE VALUE OF B DEPENDS ON THE TEMPERATURE 

8=0.0000000126 

FAC1=(PL+2. )*(2. *3. 14 16*0MEG/60 . ) 

FAC2 = B* *PL 

BETA=(FAC1/FAC2)**( 1 . /PL) 

DDX= 0.0025 

WRITE ( 6 , 567 ) B . PL , DDX , 0 , OMEG . HO 
567 FORMAT (5X.6E1 1 .4//) 

EMU= 1 . / B 
POLO= 1010000.0 
PINI =POLD 
PNEW=POLD 
PP= 0.0 
FORCE = 0.0 
INO= 0 
X = X 1 

WR I T E ( 6 , 1 00 ) XI .PP.PNEW, IND 

c MAIN LOOP, DIRECTS THE STEP-BY-STEP ADVANCEMENT OF THE SOLUTION 

DO 2 U= 1 , 800 
VIN=0/(2. *H4 ) 

POLD=PNEW 

q SECONDARY LOOP, COMPUTES THE REQUIRED DPDX VALUES AT THE INTER - 


PREOOO 1 0 

PRE00020 

PRE00030 

PRE00040 

PRE00050 

PRE00060 

PRE00070 

PRE00080 

PRE00090 

PRE00100 

PREOOI lO 

PRE00120 

PREOOI 30 

PRE00140 

PRE00150 

PRE00160 

P R E 00 170 

PREOO 1 80 

PRE00190 

PRE00200 

PRE002 10 

PRE00220 

PRE00230 

PRE00240 

PRE00250 

PRE00260 

PRE00270 

PRE00280 

PREO0290 

PRE00300 

PRE00310 

PRE00320 

PRE00330 

PRE00340 

PRE00350 

PRE00360 

PRE00370 

PRE00380 

PRE00390 

PRE00400 

PRE00410 

PRE0O420 

PRE00430 

PRE00440 

PRE00450 

PRE00460 

PR E 004 70 

PRE00480 

PRE00490 

PR E 00500 

PRE00510 

PRE00520 

PRE00530 

PRE00540 

PRE00550 
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MEDIATE GRID POINTS CONTAINED WITHIN ONE MAIN X-STEP. THE 

CALCULATED VALUES ARE STORED IN ARRAYS. 

DO 1 I =1,3 

H( I )=HO+R-SQRT(R*R-X*X) 

SL( I ) =ABS( (X4-X)/(X4+X1 ) ) 

VRX ( I )=2. *3. 1 4 1 6 *OMEG* SORT (R*R-X*X)/ 60 . 

VRXP( I )=(VRX(I )-VIN)*SL(I ) + *XN + VIN 
A ( I ) = ( PL + 2 . )*(VRXP(I)*H(I)-Q/2. )/((PL+1. ) ♦H( I ) * * ( PL + 2 . )) 
C1(I)=((PL+1 . ) * * ( 1./PL))*EMU 

DPDX ( I ) =C 1 ( I ) * ( ( ABS( A( I )))* + (( 1 . /PL) - 1 . ) ) * ( A( I ) ) 

WR I TE ( 6 , 99 ) X 

WRITE (6, 101 ) H( I ) , SL( I ) , VRX ( I ) , VRXP( I ) , A ( I ) , DPDX ( I ) 
X=X-DDX/2. 

1 ‘ CONTINUE 

END OF INNER LOOP. 

X=X+DDX/2. 

RUNGE-KUTTA FORMULA FOR THE CALCULATION OF THE PRESSURE. 


C 


99 

100 

101 

2 

C 

102 


PNEW=POLD - DDX* ( DPDX ( 1 ) + 4 . *DPDX( 2 )+DPDX( 3 ) )/6 . 
PDN 3 (PNEW-PINI ) /BETA 

IF(X. LT. ’XI +DDX . AND . X . GT . -X 1 -DDX ) PMAX= PNEW 
SHE AR= H( 2 ) *DPDX ( 2 ) 

FRICOE= SHEAR/PNEW 
FORCE = FORCE + PNEW*DDX 

IF(PNEW.LE. 1010000. 0. AND. J.GE. 25) GO TO 102 
WRITE (6 , lOO) X , DPDX ( 2 ) , PNEW , PDN , FRICOE , J 
FORMAT ( /20X , E 1 5 . 5/ ) 

FORMAT ( 1X.5E12.5.3X. 13) 

FORMAT ( 1X.6E 1 1 . 4) 

CONTINUE 

END OF MAIN LOOP. 

CONTINUE 

TORQUE = F ORCE * ABS (XI) 

FD= FORCE / ( BETA *R ) 

WRITE ( 6 , 99 ) PMAX, FORCE, TORQUE. FD 
STOP 
END 


PRE00560 
PRE00570 
PRE00580 
PRE0O590 
PRE00600 
PRE006 10 
PRE00620 
PRE00630 
PRE00640 
PRE00650 
PRE00660 
PRE00670 
PRE00680 
PRE00690 
PRE00700 
PRE00710 
PRE00720 
PRE00730 
PRE00740 
PRE00750 
PRE00760 
PRE00770 
PRE00780 
PRE00790 
PRE00800 
PRE008 1 0 
PRE00820 
PRE00830 
PRE00840 
PRE00850 
PRE00860 
PRE00870 
PRE00880 
PRE00890 
PRE00900 
PRE00910 
PRE00920 
PRE00930 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 


5.6.- Program LUDECT. The Solution of Systems of Linear 
Algebraic Equations with Tridiagonal Matrices. 


- -THIS PROGRAM SOLVES THE SYSTEM OF LINEAR ALGEBRAIC EQUATIONS 
M X = BB BY L-U DECOMPOSITION . HERE A ( N ) IS THE VECTOR 
FORMED BY THE DIAGONAL ELEMENTS OF M . B ( N ) AND C(N) ARE. 
RESPECTIVELY THE VECTORS FORMED BY THE ELEMENTS ALONG THE 
LOWER AND UPPER DIAGONALS OF M . EL ( N ) IS THE VECTOR FORMED 
BY THE ELEMENTS ON THE LOWER TRIANGULAR MATRIX AFTER L-U 
DECOMPOSITION . UP(N) IS THE VECTOR OF THE UPPER TRIANGULAR 
ELEMENTS . X(N) IS THE SOLUTION VECTOR AND N IS THE SIZE 
OF THE ORIGINAL SYSTEM . 

--THE PROGRAM IS DESIGNED SPECIFICALLY TO DEAL WITH TRIDIAGONAL 
MATRICES AND THE SPECIFIC EXAMPLE HERE IS FOR THE CASE WHEN 
M HAS 2 ' S ALONG THE MAIN DIAGONAL AND - 1'S ALONG THE 

NEIGHBORING DIAGONALS . 

--THE VECTOR ON THE RHS BB ( N ) IN THIS CASE HAS ALL COMPONENTS 
EQUAL TO ZERO EXCEPT THE FIRST . 

REAL A( 10) .B( 10) .C( lO) ,BB( 10) .EL( 10) .UP( 10) ,D( lO) 

REAL Y ( 10) ,X( 10) 

N- 10 

DO 5 1 = 1 ,N 
A ( I ) = 2. 

B( I ) = - 1 - 
C( I ) = - 1 . 

BB ( I ) -0 . 

CONTINUE 
BB ( 1 ) = 1 . 

D(l)=A(1) 

UP( 1 )=C( 1 ) 

DO 10 1=2. N 
I M 1 = I - 1 

E L ( I ) = B(I)/D(IM1) 

D( I ) = A ( I ) - EL ( I )*UP( IM1 ) 

UP ( I ) =C( I ) 

CONTINUE 


Y ( 1 ) =BB ( 1 ) 

DO 20 1=2, N 
I M 1 = I - 1 

Y ( I ) = BB ( I ) - EL( I )*Y( IM1 ) 

CONTINUE 

X(N)=Y(N)/D(N) 

DO 30 I K = 2 , N 
I -N+ 1 - IK 
IP1=I+1 

X (I )= ( 1/D( I ) )*( Y( I ) - UP ( I ) *X( IP 1 ) ) 
CONTINUE 

WR I TE ( 6 , 50 ) (X(I ) . 1=1 ,N) 

FORMAT ( 5X , F 1 5 . 8 ) 

STOP 

END 


LUDOOO 1 0 

LUD00020 

LUD00030 

LUD00040 

LUD00050 

LUD00060 

LUD00070 

LUD00080 

LUD00090 

LUD00100 

LUD001 10 

LUD00120 

LUD00130 

LUD00140 

LUD00150 

LUD00160 

LUD00170 

LUD00180 

LUD00190 

LUD00200 

LUD002 10 

LUD00220 

LUD00230 

LUD00240 

LUD00250 

LUD00260 

LUD00270 

LUD00280 

LUD00290 

LUD00300 

LUD00310 

LUD00320 

LUD00330 

LUD00340 

LUD00350 

LUD00360 

LUD00370 

LUD00380 

LUD00390 

LUD00400 

LUD00410 

LUD00420 

LUD00430 

LUD00440 

LUD00450 

LUD00460 

LUD00470 

LUD00480 

LUD00490 

LUD00500 
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5.7 


Program PS 


The Calculation of Heat Flow Including 


Convection. Patankar-Spalding Method . 


c 

c- 

c- 

c- 

c 

c- 

c- 

c 

c 

c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c- 
c - 
c- 
c- 
c- 
c- 
c- 
c 


COMMON/ I ND/R , HO , OMEG , Q , EMU . PL , TCL , TC5 . RHO , CP 

COMMON/ROLL/TROLL . HTC , TCR . ALFR 

COMMON/ ZONE /X ,X1 , X2,X3,X4.DX,DT ,XN, I VV 

COMMON/NODES/NX. NY, III 

COMMON/ COE F/ B 1 , B2 , B3 . B4 , B5 

COMMON/ V E LO/VE LX ( 12,6) , VELY( 12,6) 

COMMON/TEMP/TA ( 12,6), T( 12. 6), CR( 12, 6). RE (12. 6) 

DIMENSION TOLD( 12,6) 

THIS PROGRAM DIRECTS THE CALCULATION OF FLOW AND TEMPERATURE 

--IN LUBRICATION TYPE FLUID FLOW CONFIGURATIONS. 

TOLD IS THE INITIAL (GUESSED) TEMPERATURE FIELD. 

BOTH . NEWTONIAN AND NON-NEWTONIAN POWER LAW FLUIDS CAN 

---BE CONSIDERED . 

DATA TOLD/7 2 *760. 0/ 

DATA TOLD/ 7 2*660 . 0/ 

FOLLOWING ARE THE GEOMETRICAL AND MATERIALS DATA REQUIRED 

- R = ROLL RADIUS(CM) , HO = MI N I MUM GAP ( CM ) , OMEG = RPM OF ROLL 
-Q=VOLUMETRIC FLOW RATE PER UNIT WIDTH(CM2/S) 

- PLL = POWER LAW EXPONENT FOR LI QU ID , PLS=POWER LAW EXPONENT FOR SOL 

- EMUL = VI SCOSI TY OF THE LIQUID(G/CM S) 

-B = F LU I D I TY OF THE SOLID (UNITS DEPEND ON PLS ) 

-TCL=THERMAL CONDUCTIVITY OF THE L I QU I D ( CAL/CM S K) 

-TCS=THERMAL CONDUCTIVITY OF THE SOLID(CAL/CM S K) 

-RHO=DENSITY (G/CM3 ).CP=SPECIFIC HEAT ( CAL/G K) 

-TROLL=TEMPERATURE OF THE ROLL (K) 

-HTC=HE AT TRANSFER COEFFICIENT TO THE ROLL ( CAL/CM2 S K) 
-TCR=THERMAL CONDUCTIVITY OF THE ROLL (CAL/CM S K) 

- ALFR = THERMAL DIFFUSIVITY OF THE ROLL (CM2/S) 

-XI =LOCAT I ON OF NEUTRAL AND LIFT-OFF POINTS (CM) 

-X2*L0CAT I ON OF THE END OF SOLIDIFICATION (CM) 

-X3=L0CAT I ON OF THE BEGINING OF SOLIDIFICATION (CM) 

-X4=L0CAT I ON OF THE ENTRANCE TO THE ROLL GAP (CM) 

-NX = NUMBER OF GRIDS ALONG X, NY = NUMBER OF GRIDS ALONG Y 
-DX = GR ID SPACING ALONG X ( CM ) , DY =GR I D SPACING ALONG Y (CM) 
-XN=EXPONENT TO DESCRIBE THE SLIP IN VELOCITY AT THE SPLAT/ROLL 
-DT =F I CT I C I OUS TIME STEP REQUIRED IN THE HEAT FLOW EQUATION (S) 

R= 5.0 
HO = 0.005 
OMEG= 160.0 


PS 100010 
PS 100020 
PS 100030 
PS 100040 
PS 100050 
PS 100060 
PS 100070 
PS 100080 
PS 100090 
PS 100100 
PS1001 10 
PS 100 1 20 
PS 100 1 30 
PS 100 1 40 
PS 1 00 1 50 
PS 100160 
PS 100170 
PS 100 1 80 
PS 100190 
PS 100200 
PS 1002 1 0 
IDPS 100220 
PS 100230 
PS 100240 
PS 100250 
PS 100260 
PS 100270 
PS 100280 
PS 100290 
PS 100300 
PS 1003 10 
PS 100320 
PS 100330 
PS 100340 
PS 100350 
PS 100360 
PS 100370 
PS 100380 
PS 100390 
PS 100400 
PS 1004 10 
PS 100420 
PS 100430 


0 « 1.0 
PLS = 4.5 
B= 0.0000000126 
PLL= 1.0 
EMUL= 0.01 
TCL= O. 15 
TCS = 0.50 
RHO- 2.7 
CP= 0.25 
TROLL= 25.0 
HTC = 1.0 
TCR= 0. 1 


PS 100440 
PS 100450 
PS 100460 
PS 100470 
PS 100480 
PS 100490 
PS 100500 
PS 1005 10 
PS 100520 
PS 100530 
PS 100540 
PS 100550 
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ALFR= 0- lO 


C 

C ASSUMED VALUES OF X1.X2.X3 AND X4 

C 

COEO= Q*60./(3. 1416+OMEG) 

COE 1 = HO + R 

C0E2= SORT ( COE 1**2. - COEO ) 

C0E3= (COE 1 + C0E2 )/2 . 

X 1 = SORT ( ABS ( R * R -C0E3 + C0E3 ) ) 


C X2 = -0.57 

X2 S XI 
X3 S XI 
X4 * -0.90 
C 

C GRID PARAMETERS 

C 

NX = 12 
NY = 6 

C DXL = 0.011 

DXL = - ( X4 -X 1 ) /FLOAT ( NX - 1 ) 

C DXL = -(X4-X3)/F LOAT ( NX ) 

C DXS= 0.0633333 

DXS“ - (X2~X 1 )/FLOAT ( NX ) 

CC STEP SIZE FOR THE CALCULATION OF THE MAXIMUM LOAD QN THE ROLLS 

CC DXS= - (X4+X1 )/ ( FLOAT(NX) ) 

DT = 1.0 
C 

C EXPONENTS OF THE SLIP COEFFICIENT. 

C 

XNL= 0.5 
XNS* 0.0 
C 

C HERE ANY OF X = X4 OR X = X2 MUST BE CHOSEN ACCORDING 

C TO THE DESIRED CALCULAT I ON . I F CALCULATIONS ARE DESI- 

C -RED IN THE LIQUID REGION CHOSE THEN X = X4 . SELECT X = X2 

C --FOR COMPUTATIONS IN THE (FULLY) SOLID REGION 

C 

X= X4 
C X = X2 

C 

C THE NEXT LOOP CREATES THE INITIAL TEMPERATURE FIELD TO BE 

C SENT TO THE ROUTINE THAT SOLVES THE HEAT FLOW PROBLEM 

C 

DO 1 I = 1. NX 

DO 2 d = 1 ,NY 
T A ( I , J ) = TOLD ( I , J) 

2 CONTINUE 

1 CONTINUE 

C 

c THE FOLLOWING CONDITION COLLECTS ADDITIONAL VALUES WHICH ARE 

C FIXED AS SOON AS X HAS BEEN SELECTED. 

C 

I F ( X . GE . X2 ) GO TO 10 
9 CONTINUE 

DX= DXL 


PS 100560 
PS 100570 
PS 100580 
PS 100590 
PS 100600 
PS 1006 10 
PS 100620 
PS 100630 
PS 100640 
PS 100650 
PS 100660 
PS 100670 
PS 100680 
PS 1 00690 
PS 100700 
PS 1007 10 
PS 100720 
PS 100730 
PS 100740 
PS 100750 
PS 100760 
PS 100770 
PS 100780 
PS 100790 
PS 100800 
PS 1008 10 
PS 100820 
PS 100830 
PS 100840 
PS 100850 
PS 100860 
PS 100870 
PS 100880 
PS 100890 
PS 100900 
PS 1 009 1 0 
PS 100920 
PS 100930 
PS 100940 
PS 100950 
PS 100960 
PS 100970 
PS 100980 
PS 100990 
PS101000 
PS101010 
PS 101020 
PS 1 01030 
PS 10 1040 
PS 101050 
PS 101060 
PS 1 01070 
PS 101080 
PS 101090 
PS 101 100 
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10 


20 

c 


XN = XNL 
EMU* EMUL 
PL* PLL 
TC = TCL 
GO TO 20 
CONTINUE 
DX = DX5 
XN= XNS 
EMU* ( 1 ./B) 
PL* PLS 
TC* TCS 
CONTINUE 


WRITE (6, 99) 
WR I TE ( 6 , 99 ) 
WR I TE ( 6 , 99 ) 
WR I TE ( 6 , 98 ) 


R , HO , OMEG , 0 , EMUL , PLL, 8 , PLS. TCL, TCS, RHO , CP 
TROLL ,HTC,TCR,ALFR 
X.X1 . X2 ,X3,X4 , DX , DT , XN 
NX. NY 


THE ROUTINES FOR THE CALCULATION OF VELOCITY AND TEMPERATURE 

FIELDS ARE CALLED NOW. DEPENDING ON THE REGION OF CALCULATION 

SELECT EITHER X=X4+DX OR X=X2+DX . 

X* X4 

CALL VELOZ 

GO TO 100 
91 CONTINUE 


X= X4+DX 


DO 92 1=1 .NX 

DO 92 U= 1 , NY 
T( I , J)=TA< I , J) 

2 CONTINUE 

DO 93 I I 1 = 1 .800 

CALL TOMA 

RESMAX =0 . 0 
DO 94 IV* 1 .NX 

DO 95 UV= 1 . NY 

IF(RESMAX.LT.RE( IV, UV) ) RESMAX*RE ( IV. JV) 
RE ( IV. JV)=0.0 

5 CONTINUE 

'4 CONTINUE 

IF (RE5MAX.LT. 0.0001) GO TO 97 

3 CONTINUE 

7 CONTINUE 

WR I T E ( 6 . 99 ) ( ( T ( I , 0 ) , J* 1 , NY ) . I = 1 . NX ) 

WR I TE ( 6 , 66 ) RESMAX, III 
66 FORMAT ( // , 25X , E 1 2 . 5 . 3X , 1 3 ) 

WRITE (6, 99) R. HO, OMEG. 0, EMUL. PLL .B. PLS. TCL. TCS. RHO. CP 

98 FORMAT ( /20X .21 10/// ) 

99 FORMAT ( 1X.6E 12. 5) 

100 CONTINUE 

STOP 

END 


PS101 1 10 
PS101 120 
PS101 130 
PS101 140 
PS101 150 
PS101 160 
PS 101 170 
PS101 180 
PS 10 1 .1 90 
PS 101 200 
PS 1012 10 
PS 101220 
PS 1 0 1 230 
PS 101240 
PS 101250 
PS 1 0 1 260 
PS 101 270 
PS 101280 
PS 10 1 290 
PS 10 1 300 
PS 1013 10 
PS 101320 
PS 101330 
PS 101340 
PS 10 1 350 
PS 10 1 360 
PS 101370 
PS 10 1 380 
PS 1 0 1 390 
PS 10 1 400 
PS 1014 10 
PS 101420 
PS 101430 
PS 10 1 440 
PS 101450 
PS 10 1 460 
PS 10 1 470 
PS 1 0 1480 
PS 101490 
PS 101500 
PS 10 1 5 10 
PS 101520 
PS lO 1530 
PS 101540 
PS 10 1 550 
PS 10 1 560 
PS 101570 
PS101580 
PS 101 590 
PS 10 1 600 
PS 10 1 6 10 
PS 10 1620 
PS 1 0 1630 
PS 101640 
PS 101650 
PS 101660 
PS 1 0 1670 
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SUBROUTINE VELOZ 

THIS ROUTINE COMPUTES THE VELOCITY FIELD . 

COMMON/ I ND/R . HO . OMEG . 0 , EMU . PL . TCL . TCS . RHO , CP 
COMMON/ZONE/X ,X1 ,X2,X3.X4,DX,DT , XN , I VV 
COMMON/NODES/NX , NY 

COMMON/ VE LO/VE LX ( 12,6),VELY(12,6),Z( 12, 6), ST ( 12,6) 
COMMON/COE F/B 1 ,B2,B3,B4,B5 

DIMENSION VE LOX ( 24 , 12).VELOYf24. 12),ZINT(24. 12),STRE(24, 12) 
DIMENSION PC 26 ) 


C 

C- CALCULATION OF THE VELOCITY FIELD 

c THE VALUES OF THE FUNCTIONS THAT ONLY DEPEND ON X ARE COMPUTED 

C FIRST, WITH THE MAIN DO LOOP . THE VELOCITIES ARE CALCULATED SUB- 

C SEOUENTLY FOR FIXED X AND VARYING Y . 

C MAIN 00 LOOP 


PBR Y = 1010000.0 

DX= DX/2. 

NNX= NX * 2 - 1 

DO 1 1=1, NNX 

HI = H(R.HO.X) 

B 1 = VRXP(OMEG, R . X ,0. HO , X2 . X4.XN.X1 ) 

B2= A(R,X, OMEG , HO ,PL,Q,X2,X4, XN , X 1 ) 

B3= COE F I ( R , X , OMEG. HO . P L . 0 . X2 . X4 . XN , X 1 ) 

B4= - DADX(R,X,O t OMEG . HO ,PL,X2,X4,XN,X1)/(PL + 2. ) 

B5= VRYP ( OMEG , R , X , 0 , HO , X2 . X4 , XN, X 1 ) 

SS= 2 . *EMU*B2*HI 

DPDX = ( ( ( PL + 1 . )**(1./PL))*EMU)*((ABS(B2))**((1./PL)-1. ))*B2 
P( I ) = PBRY + DPDX *DX 

C WR I T E ( 6 . 3 ) B1 , B2 . B3 . B4 . B5 . SS , DPDX 

3 FORMAT ( 1 X , 7E 1 1 . 4 ) 

C WR I TE ( 6 , 33 ) X , DPDX . P ( I ) . HI 

33 FORMAT (5X.4E15.7 ) 

34 FORMAT (2X.5E 12.5) 

SECONDARY DO LOOP. THE VELOCITIES ARE COMPUTED FOR THE NY NODES 
LYING ALONG THE I-TH X-STEP. 


Y=0 . O 

NNY = NY *2 -1 

DO 10 d = 1 ,NNY 
H I =H ( R , HO , X ) 

D Y =H I / ( 2 . *( FLOAT ( NY -1 ) ) ) 

----CALCULATION OF VELOCITY COMPONENTS 

VX= B1 + B2* ( Y * + ( PL+ 1 . O ) - HI + * ( PL+ 1 . O ) ) 

V Y = B3* Y + B4 * Y* * ( PL + 2 . ) 

CALCULATION OF INTENSITY OF THE RATE OF DEFORMATION 

ZZ=( ABS(B2*(PL+1 ,)*Y**PL))**2. 

--CALCULATION OF THE STREAM FUNCTION 

PSI= ( ( B 1 - B2*HI +*(PL+ 1 .0) )*(HI -Y ) + 

1 (B2*(HI**(PL + 2. ) - Y**(PL + 2. ) ) )/<PL + 2. ) )*2.0/0 

IF(PSI . LE .0.0001 ) PSI=0.0 
WR I TE ( 6 , 34 ) B 1 , VX , VY , ZZ , PS I 
-THE CALCULATED VALUES ARE STORED IN ARRAYS 


PS 2000 1 0 

PS200020 

PS200030 

PS200040 

PS 200050 

PS 200060 

PS200070 

PS200080 

PS 200090 

PS 200 1 00 

PS2001 10 

PS200120 

PS200130 

PS200140 

PS200150 

PS200160 

PS200170 

PS200180 

PS200190 

PS200200 

PS2002 10 

PS200220 

PS 2002 30 

PS200240 

PS200250 

PS200260 

PS 20027 0 

PS200280 

PS200290 

PS200300 

PS200310 

PS200320 

PS200330 

PS200340 

PS200350 

PS200360 

PS200370 

PS200380 

PS200390 

PS 200400 

PS2004 10 

PS 2004 20 

PS200430 

PS200440 

PS200450 

PS200460 

PS200470 

PS200480 

PS200490 

PS200500 

PS200510 

PS200520 

PS200530 

PS200540 

PS200550 
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c 


10 
C — 


1 


VE LOX ( I , J ) = VX 
VE LOY ( I . J)=VY 
Z I NT ( I , J ) =ZZ 
STRE ( I , J ) = PS I 


Y = Y +DY 
CONTINUE 
END OF INNER LOOP 
X=X+DX 
CONTINUE 


OUTER LOOP CONCLUDED 


DO 150 1=2 .NNX. 2 

DO 150 J=1 , NNY , 2 

K=I/2 

L = (U+1 )/2 

VELX(K,L)=VE LOX( I , J) 
150 CONTINUE 


DO 160 1=3, NNX , 2 
DO 160 d = 2 , NNY , 2 
K=( 1+1 )/2 
L = J/2 

VELY ( K , L ) = VE LOY ( I , J) 
160 CONTINUE 
GO TO 180 
174 CONTINUE 

WRITING OF THE RESULTS 


WR I TE ( 6 , 177) 

177 FORMAT ( /25X , 10HVELOC I T Y -X ) 

WRITE (6, 176) ( ( VELX ( I , J ) , d = 1 , NY), 1=1, NX -1) 

175 CONTINUE 

C GO TO 180 

1755 CONTINUE 

C WR I TE ( 6 , 178 ) 

178 FORMAT ( / 25X , 1 OHVE LOC I TY - Y ) 

WR I TE ( 6 , 1766 ) ( ( VELY ( I , U ) . U= 1 , NY - 1 ) , I = 2 , NX- 1 ) 

GO TO 180 
1788 CONTINUE 

C WR I TE ( 6 , 179) 

179 FORMAT(/ 1 5X , 32HI NT ENS I TY OF RATE OF DEFORMATION) 
WRITE (6, 176) ( ( Z( I , d ) , J=1,NNY) ,1=1 .NNX) 

C WR ITE(6, 181) 

181 F0RMAT(/25X, 15HSTREAM FUNCTION) 

WRITE (6, 176) ((ST (I ,J),J=1 , NNY ) . I = 1 , NNX ) 

180 CONTINUE 

176 FORMAT ( 1X.6E12.5) 

1766 FORMAT ( 1X.5E12.5) 

C 1 80 CONTINUE 

RETURN 
END 


PS200560 
PS 2005 70 
PS200580 
PS200590 
PS 200600 
PS 2006 10 
PS200620 
PS 2006 30 
PS200640 
PS 2006 50 
PS200660 
PS200670 
PS20O680 
PS200690 
PS 200700 
PS 2007 10 
PS20O720 
PS200730 
PS200740 
PS200750 
PS200760 
PS200770 
PS200780 
PS200790 
PS 200800 
PS200810 
PS200820 
PS200830 
PS200840 
PS200850 
PS200860 
PS200870 
PS200880 
PS 200890 
PS 200900 
PS 2009 1 0 
PS200920 
PS200930 
PS2OO940 
PS200950 
PS200960 
PS200970 
PS200980 
PS200990 
PS 20 1000 
PS201010 
PS201020 
PS201030 
PS201040 
PS201050 
PS201060 
PS201070 
PS201080 
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the following functions are required for the computation 

--OF THE FLOW . 

FUNCTION H(R,HO,X) 

H=HO+R-SQRT(R*R-X*X) 

H=HO( 1 . + 0 . 05* ( X/X4 ) ) 

RETURN 

ENO 

X-COMPONENT OF THE TANGENTIAL VELOCITY OF THE ROLLS(CM/S) 

FUNCTION VRX ( OMEG , R , X ) 

VRX*2 . *3 . 1 4 1 6 * OMEG* SORT (R+R-X*X)/ 60 . 

RETURN 

END 

Y-COMPONENT OF THE TANGENTIAL VELOCITY OF THE ROLLS(CM/S) 

FUNCTION VRY ( OMEG , R , X ) 

VRY = 2 . * 3 . 14 1 6*0MEG+X/ 60 . O 

RETURN 

END 

DERIVATIVE OF H (DIMENSIONLESS) 

FUNCTION DHDX(R.X) 

DHDX = X/SQRT (R*R-X+X) 

RETURN 

END 

DERIVATIVE OF VRX (1/S) 

FUNCTION DVRXDX ( OMEG , R , X ) 

01=2. *3. 14 1 6 *OMEG/ 60 . 

DVRXDX = - C 1 * ( X/ SORT (R*R-X*X) ) 

RETURN 

END 

FLOW COEFFICIENT (UNITS DEPEND ON THE MATERIAL) 

FUNCTION A(R,X,OMEG, HO ,PL,0.X2,X4, XN , X 1 ) 

VR=VRXP( OMEG , R , X , 0 , HO , X2 , X4 , XN , X 1 ) 

HI = H(R.HO.X) 

C 1 = ( PL + 2 . )*(VR*HI-0/2. )/(PL+1 . ) 

C2 = HI * * ( PL + 2 . ) 

A= C1/C2 

RETURN 

END 

DERIVATIVE OF THE FLOW COEFFICIENT (VARIABLE UNITS) 

FUNCTION DADX(R,X,0. OMEG , HO , PL , X2 ,X4 . XN.X1 ) 

VR= VRXP ( OMEG , R , X , Q , HO , X2 ,X4 . XN.X1 ) 

VRP=DVRXPP ( OMEG , R , X . 0 , HO ,X2,X4,XN.X1 ) 

HI = H( R , HO , X ) 


PS 3000 1 O 

PS300020 

PS300030 

PS 300040 

PS 300050 

PS 300060 

PS300070 

PS 300080 

PS 300090 

P S 300 1 00 

PS 300 1 10 

PS 300 1 20 

PS300130 

PS300140 

PS 300 1 50 

PS 300 1 60 

PS 300 1 70 

PS 300 1 80 

PS300190 

PS 300200 

PS 3002 lO 

PS300220 

PS300230 

PS300240 

PS300250 

PS300260 

PS300270 

PS300280 

PS300290 

PS300300 

PS300310 

PS300320 

PS300330 

PS300340 

PS300350 

PS300360 

PS300370 

PS300380 

PS300390 

PS300400 

PS 3004 10 

PS 3004 20 

PS300430 

PS300440 

PS300450 

PS 300460 

PS300470 

PS300480 

PS300490 

PS 300500 

PS300510 

PS300520 

PS300530 

PS300540 

PS300550 


c 

c- 

c 


c 

c- 

c- 

c 


c 

c- 

c- 

c 


c 

c- 

c 


c 

c- 

c 


c 

c- 

c 


C 1 = ( PL + 2 . )/ ( ( PL + 1 . ) *HI + * ( PL + 2 . ) ) 

C2=HI*VRP+VR*DHDX(R,X) 

C3= ( PL + 2 . ) * ( VR -0/ ( 2 . *HI ) )*DHDX(R,X) 

DADX=C1*(C2-C3) 

RETURN 

END 

-FLOW COEFFICENT FOR THE Y-COMPONENT OF VELOCITY (VARIABLE UNITS) 

FUNCTION COEFI(R,X, OMEG , HO ,PL,0,X2,X4, XN . X 1 ) 

AA= A ( R , X , OMEG , HO ,PL,Q,X2,X4,XN,X1 ) 

HI* H( R , HO , X ) 

C 1 = ( (PL+1 . )*AA*HI**PL) *DHDX ( R , X ) 

C2 = DADX(R,X,0, OMEG , HO , PL t X2 , X4 , XN , X 1 )*HI**(PL+1 . ) 

C3 - DVRXPP(OMEG, R ,X , 0, HO ,X2,X4 , XN , X 1 ) 

C0EFI=C1+C2-C3 

RETURN 

END 


PS300560 
PS 3005 70 
PS300580 
PS300590 
PS 300600 
PS 3006 1 0 
PS300620 
PS300630 
PS300640 
PS300650 
PS300660 
PS300670 
PS300680 
PS 3006 90 
PS300700 
PS 3007 10 
PS300720 
PS300730 . 


PS300740 

-FUNCTION GIVING THE SURFACE X-VELOCITY OF THE SPLAT UNDER SL I PP I NGPS300750 
-CONDITIONS IN TERMS OF THE SLIP COEFFICIENT AND EXPONENT ( CM/S ) PS300760 

PS 3007 70 

FUNCTION VRXP ( OMEG , R . X t Q , HO ,X2.X4,XN.X1 ) PS300780 

SL = S(X,X2,X4 t X1 ) PS 300790 

VR= VRX ( OMEG . R . X ) PS300800 

H4 = H( R , HO , X4 ) PS3008 10 

VI N=0/ ( 2 . ♦H4 ) PS 3008 20 

VRXP=(VR-VIN)*SL**XN + VIN PS300830 

RETURN PS 3008 40 

END PS300850 

PS300860 

-FUNCTION GIVING THE SLIPPED Y-VELOCITY IN THE SURFACE OF THE SPLATPS300870 


•(CM/S) 

FUNCTION VRYP(OMEG,R , X, O.HO.X2 ,X4 , XN , X 1 ) 

SL = S ( X , X2 , X4 , X 1 ) 

VR= VR Y ( OMEG , R , X ) 

H4=H( R , HO , X4 ) 

VRYP = VR* SL * *XN 

RETURN 

END 

-SLIP COEFFICENT (DIMENSIONLESS) 

FUNCTION S ( X . X2 . X4 , X 1 ) 

S = ABS ( (X4-X- .001 )/(X4-X2- .001 ) ) 

S = ABS ( (X4-X - .001 )/(X4+X 1 -0.001 ) ) 

RETURN 

END 

-DERIVATIVE OF THE SLIP COEFFICENT (1/CM) 

FUNCTION DSDX ( X2 , X4 , X 1 ) 

DSDX= - 1 ./(-X4-X2) 

DSDX* - 1 ./(X4 + X1 ) 

RETURN 

END 

-FUNCTION GIVING THE X-DERIVATIVE OF THE SLIPPED VELOCITY (1/S) 

FUNCT I ON DVRXPP ( OMEG , R , X . 0 , HO , X2 . X4 . XN , X 1 ) 

SL = S(X,X2,X4,X1 ) 

VIN=0/(2. *H( R , HO , X4 ) ) 

VR = VRX ( OMEG , R , X ) 

OVRXPP = ( ( VR -VIN) *XN ) *(SL**(XN- 1 . ) ) *DSDX ( X2 . X4 , X 1 ) 

1 +(SL**XN)*( DVRXDX( OMEG , R , X ) ) 

RETURN 

END 


PS300880 
PS300890 
PS 300900 
PS 3009 10 
PS300920 
PS300930 
PS300940 
PS300950 
PS300960 
PS300970 
PS300980 
PS300990 
PS 30 1000 
PS301010 
PS301020 
PS301030 
PS301040 
PS301050 
PS301060 
PS301070 
PS301080 
PS 30 1090 
PS301 100 
PS301 1 10 
PS301 120 
PS301 130 
PS301 140 
PS301 150 
PS301 160 
PS301 170 
PS301 180 
PS 30 1 190 
PS 30 1 200 
PS3012 10 
PS301220 
PS301230 
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c- 

c- 

c- 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c- 

c 

c 

c 

c 

C2 

Cl 

c 

c- 

c 

c 


SUBROUTINE TOMA 

--THIS ROUTINE USES A TWO-DIMENSIONAL VELOCITY FIELD TO 
COMPUTE THE CORRESPONDING TEMPERATURE FIELD. 

--THE ROUTINE ASSUMES NO HEAT FLOW BOTH, ACROSS THE SYMMETRY 
PLANE ( Y =0 . 0 ) AND THE LOWER HORIZONTAL PLANE (X = X3 OR XI) AS 
WELL AS PRESCRIBED TEMPERATURE AT THE UPPER HORIZONTAL BOUN- 
DARY ( X - X4 OR X=X2).AND PRESCRIBED HEAT FLUX TO THE ROLLS. 

--THE CALCULATION PROCEEDS BY SOLVING THE SYSTEM OF NX 
EQUATIONS WITH NY UNKNOWNS EACH. OBTAINED FROM F INI TE -D I FFEREN- 
CING OF THE HEAT FLOW EQUATION. THE METHOD IS ITERATIVE INASMUCH 
•AN ASSUMED TEMPERATURE FIELD IS USED TO COMPUTE AN IMPROVED 
■GUESS WHICH IN TURN IS USED TO COMPUTE AN EVEN BETTER GUESS. 
■THE PROCEDURE IS REPEATED UNTIL SATISFACTORY CONVERGENCE IS 
■REACHED . 


COMMON/ I ND/R . HO . OMEG , Q, EMU , PL . TCL . TCS , RHO , CP 
COMMON/ROLL/TROLL, HTC ,TCR . ALFR 
COMMON/ ZONE/X , X 1 ,X2,X3,X4,DX,DT,XN, IVV 
COMMON/NODES/NX , NY . I I I 

COMMON/ VELO/VELX( 12,6) .VELY( 12,6) 

COMMON/ TEMP/T A ( 12.6).T( 12.6) ,CR( 12.6) .RE ( 12,6) 
COMMON/ I NTCOE/COO .CIO. C20 , C01 2 


THE FOLLOWING ARRAYS ARE NECESSARY FOR THE SOLUTION 

THE SYMBOLS ARE THE SAME AS IN S . PATANKAR , NUMERICAL 

HEAT TRANSFER AND FLUID FLOW, HEMISPHERE .WASHINGTON. 1980. CH . 5 . 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
NXN = NX - 1 
NYN=NY - 1 
TCL= O. 15 
TCS= 0.52 
CPL= 0.25 
CPS= 0.25 


A( 12,6),B( 12,6),C( 12,6),D( 12,6),P( 12,6),Q( 12,6) 
DN( 6) ,DS(6) ,DE(6) , DW( 6 ) 

PEN(6),PES(6) , PE E ( 6 ) ,PEW(6) 

AN 1 ( 6 ) ,AS1(6) .AE1 (6) . AW 1(6) 

AN2( 6) ,AS2(6) ,AE2(6) ,AW2(6) 

COEN( 6 ) , COES( 6 ) , COEE (6) ,C0EW(6) 

TG( 12,6) 


THE FOLLOWING LOOP ONLY INITIALIZES THE TEMPERATURE FIELD 


DO 1 1=1. NX 

DO 2 J= 1 , NY 
T ( I , J)=TA( I , J) 

CONTINUE 

CONTINUE 

THE FOLLOWING LOOP CONTROLS THE ITERATIONS 

DO 3 111=1, 800 
C00 = 25 . 

C10=0. 

C20 = 0 . 0 


PS4000 10v 

PS400020 

PS400030 

PS400040 

PS 400050 

PS400060 

P S 400070 

PS 400080 

PS400090 

PS400100 

PS 400 1 10 

PS400120 

PS400130 

PS400140 

PS400150 

PS400160 

PS 400 170 

PS400180 

PS400190 

PS400200 

PS4002 10 

PS400220 

PS400230 

PS400240 

PS400250 

PS400260 

PS400270 

PS400280 

PS400290 

PS 400300 

PS400310 

PS400320 

PS400330 

PS400340 

PS400350 

PS400360 

PS400370 

PS400380 

PS400390 

PS 400400 

PS4004 10 

PS400420 

PS400430 

PS400440 

PS400450 

PS400460 

PS400470 

PS400480 

PS400490 

PS400500 

PS400510 

PS400520 

PS400530 

PS400540 

PS400550 
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c- 

c- 

c- 

c- 

c 

cc 

c 

c- 

c 


c 

15 

c 

c- 

c- 

c- 

c- 

c 


c 

c 

c 

c 


c 

c 

c 

c 


THE FOLLOWING LOOP ADVANCES THE SOLUTION IN THE X DIRECTION. 

THE COMPUTATION BEGINS AT 1=2 BECAUSE AT 1=1 THE TEMPERATURES 
ARE GIVEN AS BOUNDARY CONDITION BY THE POURING TEMPERATURE OR 
THE SOLIDUS TEMPERATURE . DEPEND I NG ON THE REGION OF COMPUTATION. 

DO 10 I - 2 , NXN 

DTI = 100. * ( 60. / ( 2 . *3 . 14 16'OMEG) ) * ( DX/SQRT ( R*R- ( X -DX ) * ( X -OX ) ) ) 

THE NEXT LOOP SETS THE INLET TEMPERATURE. 

DO 15 J= 1 , NY 
T ( 1 , J)=TA( 1 . J) 

T( 1,d)*TA( 1,J) - 4 . 0*FL0AT ( J ) + 4. 

CONTINUE 

WE BEGIN NOW THE CALCULATION OF THE COEFFICIENTS FOR THE HEAT 
•FLOW EQUATION (DISCRETIZED). 

THE NODES ALONG THE SYMMETRY LINE (Y=0.0) HAVE COEFFICIENTS 

GIVEN BY THE SYMMETRICAL BOUNDARY CONDITION. 

A ( I , 1 ) = 1 . 

B( I , 1 )=1 . 

C ( I , 1 ) =0 . 

D( I . 1 )=0. 

CALCULATION OF THE COEFFICIENTS FOR GAUSS ELIMINATION 

FOR THE NX NODES LYING ALONG THE SYMMETRY LINE(I.I) 

P ( I , 1 )=B( I . 1 ) / A ( 1,1) 

Q( I . 1 )=D( I , 1 ) / A ( 1,1) 

-WE CONTINUE NOW WITH THE SECONDARY LOOP WHICH ADVANCES THE CAL- 
CULATION OF THE COEFFICIENTS ALONG THE Y-OIRECTION. 

DO 20 J = 2 , NYN 


C-- THE NEXT CONDITION IS REQUIRED TO MAINTAIN THE VALUES OF TC 

C AND CP THROUGHOUT THE CALCULAT ION . I T ENDS IN LABEL 5. 

C 

IF (X.GE.X2) GO TO 4 
TC= TCL 
CP= CPL 
GO TO 5 

4 CONTINUE 
TC = TCS 
CP= CPS 

5 CONTINUE 

HI =H( R , HO , X ) 

DDY = HI / (FLOAT (NYN) ) 

IF(J.GE.NY) GO TO 21 
JJ=J+ 1 
JJU=J- 1 
C 

C CALCULATION OF THE COEFFICIENTS FOR THE INTERNAL GRID POINTS. 


PS 4005 60 

PS 4005 70 

PS400580 

PS400590 

PS 400600 

PS400610 

PS400620 

PS400630 

PS400640 

PS400650 

PS400660 

PS400670 

PS400680 

PS400690 

PS400700 

PS400710 

PS400720 

PS400730 

PS400740 

PS400750 

PS400760 

PS400770 

PS400780 

PS400790 

PS400800 

PS4008 10 

PS400820 

PS400830 

PS 4008 40 

PS400850 

PS400860 

PS400870 

PS400880 

PS400890 

PS400900 

PS 4009 10 

PS400920 

PS400930 

PS400940 

PS4OO950 

PS400960 

PS400970 

PS400980 

PS400990 

PS 40 1 000 

PS401010 

PS401020 

PS401030 

PS401040 

PS401050 

PS401O60 

PS401070 

PS401080 

PS401090 

PS401 100 
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c COEFFICIENTS ON THE NORTH SIDE (Y+). PS401110 

C PS401120 

DN(d)=TC*DX/(CP*DDY) PS401130 

PEN(d)=ABS(VELY( I , dd ) *DDY »RHO*CP/TC ) PS401 140 

AN 1 ( 0 ) = AMAX 1(0. , ( ( ABS ( 1.-0. 1 *PEN( d ) ) ) * *4 . ) + ( 1 . -O. 1*PEN(d) ) ) PS401 150 
AN2 ( d ) = AMAX 1 ( -RHO«DX*VELY( I ,JJ) .0. ) PS401 160 

C PS401170 

C COEFFICIENTS ON THE SOUTH SIDE (Y-). PS401180 

C PS401190 

OS(d)=TC*DX/(CP*DDY) PS401200 

PES(d)=ABS(VELY(I , JJJ ) *DDY«RHO*CP/TC ) PS401210 

AS 1 ( J ) =AMAX 1(0. . ( ( ABS( 1.-0. 1 *PES( d ) ) ) * *4 . )* ( 1 . -0. 1 *PES( J ) ) ) PS401220 
AS2 ( d ) = AMAX 1 (RHO*DX*VELY( I , ddd) ,0. ) PS401230 

C PS401240 

C COEFFICIENTS ON THE EAST SIDE (X+). PS401250 

C PS401260 

DE(d)=TC*DDY/(CP*DX) PS401270 

PEE(d)=ABS(VELX(I+1 . J ) »DX*RHO*CP/TC ) PS401280 

AE 1 ( J ) = AMAX 1(0. . ( ( ABS( 1 . -0 . 1 +PEE ( d ) ) ) * *4 . ) ♦ ( 1 . -0 . 1 *PEE ( d ) ) ) PS401290 
AE 2 ( d ) = AMAX 1 ( - RHO *DD Y * VE LX ( 1 + 1 .0) .0. ) PS401300 

C PS401310 

C COEFFICIENTS ON THE WEST SIDE (X-). PS401320 

C PS401330 

DW(J)=TC*DDY/(CP*DX) PS401340 

PEW(d)=ABS( VELX( I - 1 , d ) *DX * RHO*CP/TC ) PS401350 

AW 1 ( <J ) = AMAX 1(0. . ( ( ABS ( 1.-0. 1*PEW(d)))**4. )♦( 1.-0. 1 *PEW( d ) ) ) PS401360 
AW2 ( d ) = AMAX 1 ( RHO *DDY * VE LX ( I - 1 .d) .O ) PS401370 

C PS401380 

c CALCULATION OF GLOBAL COEFFICIENTS FOR N.S.E.AND W PS401390 

C PS401400 

COEN( d ) =DN (d)*AN1(d) + AN2 ( d ) PS401410 

C0ES(d)=DS(d)*AS1(d) + AS2(d) PS401420 

COE E (d)=DE(d)*AE1(d) + AE2( d ) PS401430 

C0EW(d)=DW(d)*AW1(d) + AW2(d) PS401440 

A0 = RHO *DX *DD Y /DT PS401450 

C PS401460 

C CALCULATION OF THE MAIN COEFFICIENTS. PS401470 

C PS401480 

A ( I ,d)=COEN(d)+COES(d)+COEE(d)+COEW(d)+AO PS401490 

B( I , d ) =COEN( d ) PS401500 

C ( I ,d)=COES(d) PS401510 

D( I , d )=COEE(d)*TA( 1 + 1 , d)+COEW(d)*T( I- 1 , d ) + AO*TA ( I . d ) PS401520 

C PS401530 

c CALCULATION OF THE COEFFICIENTS FOR THE GAUSS ELIMINATION PS401540 

c FOR THE INTERNAL NODES ( I =2 . NX - 1 ; d = 2 . NY - 1 ) . PS401550 

C PS401560 

P( I ,d)=B( I ,d)/(A( I ,d)-C(I ,d)«P( I , ddd ) ) PS401570 

Q( I ,d)=(D( I ,d)+C( I ,d)*Q( I .ddd) )/( A( I,d)-C(I.d)*P(I. ddd) ) PS40 1 580 
C PS401590 

20 CONTINUE PS401600 

C PS401610 

C END OF INNER LOOP. PS401620 

C PS401630 

21 CONTINUE PS401640 

C PS401650 
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C PS401660 

C THERMAL BOUNDARY CONDITION AT STRIP/ROLL INTERFACE. PS401670 

CC INTRODUCTION OF IDEAL COOLING. PS401680 

CC IVV=I-1 PS401690 

CC CALL IDEALC PS401700 

CC A ( I , NY ) = TC/DDY + TCR/SQRT ( 3 . 1 4 16 * ALFR *OT I ) PS401710 

CC D( I ,NY)=TCR*C012/SQRT(3. 14 16*ALFR*DTI ) PS401720 

C PS401730 

C CALCULATION OF THE EFFECTIVE HEAT TRANSFER COEFFICIENT PS401740 

C , PS401750 

VR=VRX(OMEG,R,X) PS401760 

VRP=VRXP ( OMEG . R , X , 0 . HO , X2 , X4 , XN. X 1 ) PS401770 

HTCP = HT C*(VRP/VR) PS401780 

C PS401790 

C USE OF THE HEAT TRANSFER COEFFICIENT FOR BOUNDARY NODES PS401800 

C AT THE SPLAT/ROLL INTERFACE. PS401810 

C PS401820 

A (I , NY ) = TC/DDY + HTCP PS401830 

B ( I , NY ) =0 . PS401840 

C( I .NY )=TC/DDY PS401850 

D(I.NY)* HTCP*TROLL PS401860 

C PS401870 

c AT THIS STAGE. THE VALUES OF ALL THE COEFFICIENTS FOR THE NY PS401880 

C Y-NODES LYING AT THE I-TH X-STEP.HAVE BEEN CALCULATED. PS401890 

C PS401900 

C EVALUATION OF THE N-TH COEFFICENT FOR THE GAUSS ELIMINATION. PS401910 

C COEFFICIENT MAKE EQUAL TO THE TEMPERATURE OF NODE I, NY. PS401920 

C PS401930 

Q( I ,NY) = (D(I ,NY)+C(I . NY ) *Q( I , NYN ) )/ ( A( I . NY ) -C( I , NY ) *P ( I , NYN ) ) PS401940 

T(I ,NY)=Q( I .NY) PS401950 

C PS401960 

C THE FOLLOWING RECURSIVE NODE COMPUTES TEMPERATURES FOR NODES PS401970 

C LYING ALONG THE I-TH X-STEP, TRAVELING FROM THE SPLAT/ROLL PS401980 

C INTERFACE TOWARDS THE CENTER OF THE SPLAT. PS401990 

C PS402000 

DO 30 K= 1 , NYN PS402010 

KK=NY -K PS402020 

KKK=NY+ 1 -K PS402030 

T( I , KK ) = P ( I , KK ) * T ( I . KKK ) + Q(I.KK) PS402040 

30 CONTINUE PS402050 

C PS402060 

C -THE FOLLOWING LOOP COMPUTES THE RESIDUAL BETWEEN THE PS402070 

C FRESHLY CALCULATED TEMPERATURE FIELD AND THE ONE OBTAINED PS402080 

C -FROM THE PREVIOUS ITERATION. PS402090 

C PS402100 

DO 35 KL= 1 , NY PS4021 10 

RE ( I , KL ) = ( ABS ( T A ( I , KL ) -T < I . KL ) ) )/TA ( I , KL ) PS402120 

35 CONTINUE PS402130 

C PS402140 

C THE NEXT LOOP CREATES A NEW (IMPROVED) GUESS FOR PS402150 

C THE TEMPERATURE FIELD FROM THE FRESHLY CALCULATED. PS402160 

C PS402170 

DO 40 L = 1 , NY PS402 180 

TA( I . L)=T( I . L ) PS402190 

40 CONTINUE PS402200 
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IF( I -GE .NX) GO TO 60 
10 CONTINUE 

C 

c end OF MAIN (X) LOOP . AT THIS STAGE THE ENTIRE TEMPERATURE FIE 

c IS KNOWN EXCEPT THE VALUES FOR THE NY GRID POINTS LYING ALONG 

c the NX-TH x-step. 

C 

60 CONTINUE 

C 

q ---THE NEXT LOOP USES THE BOUNDARY CONDITION OF ZERO HEAT FLOW 

c ALONG THE X-DIRECTION AND COMPLETES THE TEMPERATURE FIELD. 

C 

DO 80 LL= 1 .NY 

C T ( NX , LL ) - 660. 

c TA(NX,LL)= 660. 

T(NX,LL)=T( nxn . LL ) 

ta(nx,ll)=ta( nxn , LL ) 

RE ( NX , LL ) =RE ( NXN , LL ) 

80 CONTINUE 
GO TO 700 

81 CONTINUE 
C 

C- THE NEXT LOOP SCANS THE ENTIRE TEMPERATURE FIELD COMPARING 

C -THE FRESH VALUES WITH THOSE OBTAINED FROM THE PREVIOUS ITE- 

<* RATION AND LOOKS FOR THE LARGEST RESIDUAL. 

C 

C RESMAX =0 . O 

C DO 85 I V = 1 , NX 

C DO 90 JV=1 .NY 

C I F ( RESMAX . LT . RE ( I V . JV ) ) RESMAX-RE ( I V . JV ) 

C RE ( I V , JV ) =0 . 0 

90 CONTINUE 

85 CONTINUE 

C 

C SATISFACTION OF THE CRITERION FOR CONVERGENCE. 

C 

c I F ( RESMAX . LT . 0 . 000 1 ) GO TO 69 

03 CONTINUE 

C 

69 CONTINUE 

C69 X3=X 

C 

C WRITING OF THE CONVERGED TEMPERATURE FIELD. 

C 

WRITE (6 . 100) ( (T( I , J) . J= 1 .NY ) , 1= 1 .NX) 

C 

0 CALCULATION OF COOLING RATE FIELD. 

C 

DO 76 I K = 2 , NX 
I IK = IK - 1 

DO 75 JK - 1 .NY 

TG( IK,JK)=(TA( IK,JK)-TA( IIK.JK) )/DX 
CR( IK,JK) = ABS( (TG( IK.JK))*(VELX(IK.JK)) ) 

75 CONTINUE 

76 CONTINUE 

0 

0 WRITING THE FINAL COOLING RATE FIELD. 

WRITE (6. 100) ((CR(I,J),U=1,NY).I=1.NX) 

WR I TE ( 6 , 66 ) RESMAX, III 
66 FORMAT (//.25X.E12.5.3X, 13) 

100 FORMAT ( 1X.6E12.5) 

700 RETURN 

END 


PS402210 
PS402220 
PS402230 
LDPS402240 
PS402250 
PS402260 
PS402270 
PS402280 
PS402290 
PS402300 
PS402310 
PS402320 
PS402330 
PS402340 
PS402350 
PS402360 
PS402370 
PS402380 
PS402390 
PS402400 
PS4024 10 
PS402420 
PS402430 
PS402440 
PS402450 
PS402460 
PS402470 
PS402480 
PS402490 
PS402500 
PS402510 
PS402520 
PS402530 
PS402540 
PS402550 
PS402560 
PS402570 
PS402580 
PS402590 
PS402600 
PS402610 
PS402620 
PS402630 
PS402640 
PS402650 
PS402660 
PS402670 
PS402680 
PS402690 
PS402700 
PS4027 10 
PS402720 
PS402730 
PS402740 
PS402750 
PS402760 
PS402770 
PS4027S0 
PS402790 
PS402800 
PS402810 
PS402820 
PS402830 
PS402840 
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APPENDICES 


2A.1.- A Comment on the Relationship Between the Size of 

Microstructural Features and the Casting Parameters, 

Many rapidly solidified samples have been observed to have 
microstructures reminiscent of the well known cellular-dendritic 
structures typical of more conventionally cast samples. The 
main difference is that the size of these microstructural fea- 
tures becomes smaller the larger the rate of heat extraction 
from the sample. During the past few decades metallurgists have 
been searching for appropriate correlations capable of represent- 
ing the relationship between the microstructural features and 
the solidification parameters. A convenient measure of the effect 
of casting parameters on the material microstructure has been 
found to be the so called dendrite arm spacing (i.e. the center 
to center distance between neighboring dendrite arms) . Both 
primary and secondary dendrite arm spacings have been widely 
used . 

It is commonly believed that the dendrite arm spacing adjusts 
itself to the prevailing growth conditions by reducing the 
supercooling of the liquid lying between the dendrites to a 
minimum value. The experimental results indicate that the prod- 
uct of the temperature gradient in the liquid in front of the 
solid-liquid interface, , and the growth velocity , R , 

correlates well with the measured arm spacings. The proposed 
correlation has the following form , 
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n 


(i) 


X 



Note that , in Eqn(l) , the product G R has the units of 
cooling rate (i.e. °C/s) and in this sense it can be considered 
as an "effective" cooling rate just ahead of the solidification 
interface, i.e. 



Although the phenomenon of dendrite arm coarsening has been 
mentioned as one of the reasons for the occasional discrepancy 
found between Eqn(l) and actual measurements, this process is 
not believed to be very important during RS because of the 
small time intervals involved in the completion of the solid- 
ification. Moreover, in rapidly solidified melt spun samples, 
primary and secondary dendrite arm spacings are not easily 
recognized (Speck(1985) ) , and what one usually sees are cell- 
like structures which are very fine in the portion of the splat 
nearest to the wheel surface and which increase in size with 
increasing distance from this surface. It has also been suggested, 
(Kattamis(1981)) , that the microstructural sizes measured in such 
samples be considered analogous to the secondary dendrite arm 
spacings observed in the more conventionally cast samples. 
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The rapid solidification of eutectic alloys by coupled eutectic 
growth produces a lamellar microstructure much finer than the one 
found in conventionally cast samples. In this case the inter- 
lamellar spacing is found to be given by 



- 1/2 

B 3 ( R ) 


(3) 


Equations (1) and (3) above can be used to predict the scale of 
the microstructure for a wide variety of alloys once the solid- 
ification conditions have been defined. The required values of the 
coefficients and for several alloy systems can be 

obtained from Table(l) below. 


153 



Table (2A.1.1).- Parameters for the Calculation of Cell-Secondary 
Dendrite Arm Spacings from the Solidification Conditions. See Eqns 
(1) and (3). From Jones(1982). 


Alloy 

Cooling Rate Range 


B 1 

B 3 

n 

(wt%) 

(K/s) 

(y*m 

./(K/s)" 0 ) 

,, .3/2, 1/2. , . 

((M.m) / s ) (-) 

Sn-15Pb 

0.005 - 50 


23 

- 

0.35 

Al-4 . 5Cu 

0.00002 - 300 


41 

- 

0.39 

Al-10 . 5Si 

400 - 1.2*10 5 


47 

- 

0.33 

Cu-0.5Zr 

1 - 1*10 7 


160 

- 

0.40 

Inconel 

718 

0.1 - 100 


34 

- 

0.34 

X-40 

(Co-Cr) 

0.1 - 200 


40 

- 

0.27 

Fe-20Mn 

60 - 1400 


150 

- 

0.25 

Fe-25Ni 

0.001 - 1.7*10 6 


60 

- 

0.32 

440 S. Steel 15 - 1*10 5 


60 

- 

0.41 

Maraging 300 

Steel 0.1 - 1400 


40 

- 

0.30 

Ti- 2 to 30 
Al,Fe,Ge,Mo 

or V 0.2 - 150 


16 - 124 


0.3 - 0 
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Table (2A.1.1).- (contd.) 


Alloy 

(wt%) 

A1-A1 2 Cu 

Al-Zn 

Bi-Zn 

Cd-Zn 

Pb-Ag 

Pb-Cd 

Sn-Ag 

Sn-Cd 

Sn-Pb 

Sn-Zn 


Cooling Rate Range 
(K/s) 

not given 


B, 


(y*.m/( K/s) n ) (ymm ) 3/2 / 


10.5 - 11. 

8.0 

8.3 

5.3 

11.0 

4.5 
16.7 

8.5 

5.5 

8.3 


n 

s l/2> (-) 

8 
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3A.1,- The Governing Equations of Transport Phenomena, 
Introduction 

The governing equations of transport phenomena simply express 
well known facts about nature in mathematical symbols. These facts 
are based on physics, specifically, Newton's laws and Thermodynamic 
principles. There are two alternative ways of deriving the 
governing equations for momentum, heat and mass transfer in 
continua, namely; 

(i) differential (shell) balances (e.g. Bird et al(1960) and 
Szekely and Themelis(1971) ) . In these one uses the basic laws 
of physics to establish mass, momentum and energy balances over 
a small volume element inside the material under study. The 
partial differential equations resulting from making the volume 
of the element go to zero constitute the differential balances 
for all points inside the domain. 

(ii) The postulat ional approach (e.g. Slattery (1981) and Billington 
and Tate(1981)). Here one starts from a small set of postulates 
(also based on basic physics) and uses a mathematical result known 
as the transport theorem to arrive to the governing equations. 

These equations are general and valid for all materials and must 

be specialized by introducing constitutive equations of material 
behavior to arrive at the forms most useful for applications. 
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For the sake of completeness and since the postulational approach 
does not seem to have been widely known among metallurgists, we 
have decided to present here a brief summary of the method. After 
the derivation of the general equations for transport phenomena 
we will then discuss the boundary conditions most frequently 
found in fluid flow and heat and mass transfer problems. The 
section concludes with a comment on the special form the equations 
have for the case of the flow of thin liquid films and also the 
case of flows with negligible inertia. For details, however, the 
references given above should be consulted. 

The Transport Theorem and the Postulates of Continuum Mechanics 

The transport theorem is simply a mathematical identity which 
can be proved to be satisfied by any scalar, vector, or tensor 
valued function of time and position. We will restrict ourselves 
to the presentation of two of the most useful forms of the 
theorem. 

If § is a scalar, vector, or tensor and V is the volume 
of the material under study, the following is true. 


.(/$ 


d( J V dV)/dt = 

V V 


/ ( D J/Dt + £> 


/Dt + ©div v ) dV 


( 1 ) 


and , if mass is conserved (i.e. if div v_ - 0 ), 


d( // £dV)/dt - /(P D $/Dt) dV 
V V 


( 2 ) 


157 


The postulates of continuum mechanics are simply the summary of 
centuries of empirical observation of physical processes taking 
place in the world. They could be listed as follows; 

a) the principle of conservation of mass: " The mass of a given 

body is independent of time". In symbols 


d ( 



0 


(3) 


b) The principle of conservation of momentum: "The rate of change 
of momentum of a body is equal to the sum of the forces acting on 
it", i.e. 


J 


j) v dV)/dt 



(A) 


where T ' n is the field of contact forces. 

c) The principle of conservation of energy: "The rate of change of 
the total energy of a body is equal to the rate of work done on the 
body plus the rate of energy transmission to it", i.e. 


b 

TT * 


d( J P (U + v /2) dV 

V J A 


= J v * (T ’n) dA + f f (v*f) dV + 


+ / h dA + J P q 


p Q dV 


( 5 ) 


where the terms on the R.H.S. are , respectively, the rate of 
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work done by the contact forces, the rate of work done by the 
external forces, the rate of energy transmission from the 
enviroment to the body through its surface, and the rate of 
energy generation inside the body. 

Equations (3) -(5) are valid regardless of the size of the body, 
thus , the relations valid for integrated quantities are also 
valid for the quantities inside the integral signs. If one 
uses now the transport theorem, the following forms of the 
equations can be readily derived. 


a) conservation of mass 


f / 3 t + div( Ji v ) = 0 


( 6 ) 


b) conservation of momentum 


+ (7v)‘v = divT + yOf_ 


(7) 


and (c) conservation of energy. 


f> ?0/^t + (7 D)‘ v ■ - div q + tr( T ' 7v ) + 


+ f q 


( 8 ) 
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Since Eqns(6)-(8) apply regardless the material constitution 
of the body in question, constitutive equations for material 
behavior must now be introduced to account for the widely different 
properties of various materials. The constitutive equations 
are simply relationships between quantities in Eqns(6)-(8) which 
make the transport problem well posed. Typically, the flux of 
momentum T^ , and the flux of energy, q , are related to the 
intensity driving momentum transfer, ^7 X > and to the intensity 
driving heat flow, ^ u , respectively . 

Constitutive Equations for Material Behavior 

The constitutive equations are relationships between fluxes 
and intensities which enable us to classify materials into various 
groups based on similarities in their mechanical and thermal 
behaviors. Several rules have been laid down for the construction 
of appropriate constitutive equations. However, most of the most 
useful ones are based on the results of a great deal of empirical 
research. Here we restrict ourselves to presenting some of the 
most widely used constitutive equations. 

a) The Newtonian fluid : ,r The stress is a linear function of the 
strain rate M , i.e. 

^ = - pi + 2M(( y v) + ( yv) T )/2 (9) 

b) The non-Newtonian, power law fluid -creeping solid, 
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( 10 ) 

(10a) 
(10b) 

where B^g and b are material properties and can be 

obtained for the case of creeping solids from the compilation by 

Frost and Ashby(1982). Moreover, since D and S .. , the 

err ert 

effective deformation rate tensor and deviatoric stress tensor, 
are functions of the strain rate, in Eqn(10) is clearly 

a non linear function of the strain rate. 

For heat transfer, the most widely used constitutive equation 
is Fourier's first law. This is, 

q = - K(u) V u (11) 

The substitution of Eqns(9) and (11) into (6)-(8) leads to the 
forms of the governing equations which are the strating point 
for the calculations in this thesis. I.e. 

a) the equation of continuity 

'bf/'b* + div( jO v^ ) = 0 (12), 


pi + 2 y (( Vv) + (7v) T )/2 


where 


and 


7 - 


S eff /3 D eff 


D eff B 10 ( S eff ) 
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b) the equation of motion (Navier-Stokes) , 


f'bvj'b t + ( ^ v ) * V = - <7p +y A div(yv) + (13) 

and (c) the differential energy balance, 

‘JE/'Jt + (VE)’v = div ( K ^ u) + tr((T + pX) ’ V v) + 

+ J> Q (14) 

Equations (12)-(14) are the starting point of all of our 
calculations. However, before embarking oneself on the problem 
of solving these equations one has to look for suitable 
boundary conditions. Moreover, multicomponent systems are very 
common in practice and due account must be taken for them. So, 
in the next section we comment on the formulation of mass 
transfer problems in multicomponent systems and then we discuss 
the most commonly used boundary conditions for transport phenomena 
problems . 

Multicomponent Systems: Mass Transfer 

Although the presence of several components complicates the 
situation, relatively few additional ideas are required to 
formulate the governing equations for the case of multicomponent 
systems. Specifically, each of the component species in the 


162 


system can be regarded as a continuous medium with a variable mass 
density field. The model for the multicomponent mixture is then a 
superposition of all these continuous media. 

A new vector field is also introduced to describe the rate of 
motion of each species in the system. Such quantity is called the 
mass flux vector j . The mathematical indeterminancy resulting 
from the introduction of the mass flux vector has to be resolved 
by establishing constitutive relationships linking j with 
the corresponding intensities for mass transfer. This is 
reminiscent of the procedure followed before to transform Eqns(7) 
and (8) into (13) and (14) . The simplest possible case is that 
of binary systems. In this case the constitutive equation is 
known as Fick's first law. For the binary system, the equation of 
continuity takes the following form, 

^C/?t + ( Vcr V = div( D yc) + r c (15) 


Boundary Conditions for Momentum, Heat and Mass Transfer 

For Eqns(12)-(15) to constitute a well posed problem, boundary 
and initial conditions representing the specific systems under 
study must be added. These conditions are nothing but restrictions 
on the values of the field variables or their fluxes, imposed by 
the actual physics of the problem, at some locations on the 
computational domain. We now describe some of the most widely used 
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boundary conditions for the solution of transport phenomena 
problems . 

a) Boundary conditions for fluid flow problems. 

(i) Continuity of tangential components of velocity at phase 
interfaces. This is the well known no-slip condition (Batchelor 
(1967)). In symbols, denoting by a subindex a, b the phases 
in question, we have 


( * ' i > (16). 

a b 

(ii) Discontinuity of tangential components of velocity at phase 
interfaces . This condition, as opposed to the previous one, allows 
for some slip at the interface between phases in relative motion. 

It seems to be particularly useful when dealing with problems 
involving lines of contact (e.g. Dussan(1979) ) . In symbols, one 
possible way of expressing this, is 


( v * t ) =-K(v*_t) 

a b 


(17) . 


(iii) Continuity of stress at phase interfaces. Both the 
tangential and the normal components of the stress acting on the 
phase interface must satisfy continuity relationships which 
involve stresses due to surface tension effects ( Levich and 
Krylov(1969) ) . For an interface separating two Newtonian fluids, 


16 4 


the balance of tangential stresses is 


fc i n j > = y* D ij n j > b + ^ ^ i as) 

where = (( V v) + (V v)^ )/2 , and the usual 

summation convention of adding over repeated subscripts has been 
used . 

And the balance of normal stresses is, 



♦ <T(2^) 


( 19 ). 


b) Boundary Conditions for Heat Transfer Problems. 

(i) Temperature specified at the interface (Dirichlet condition) 
In symbols 

u - u (20) . 

a b 

(ii) Heat flux specified at the surface (Neumann condition) .Here, 

( K y u ) n = f(x,t) (21). 

a 

(iii) Heat flux at the surface specified by a heat transfer 
coefficient (Cauchy condition). In this case. 
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( K ^7 u ) • n = - h ( u - U-0 ) (22). 

a 

iv) Continuity (or discontinuity) of the heat flux at the phase 
interface. This condition is also called ideal cooling (in the 
case of continuity) or the Stefan condition (in the case of 
moving boundary problems with a discontinuity) . For ideal 
cooling , in symbols, 

( K V u ) = ( K V u ) ( 23 ) . 

a b 

When heat is released/absorbed at the phase interface due, for 
example to a change of phase (Stefan problem), this must be 
accounted for by adding a corresponding term to Eqn(23), (see 
Sec(3A.2) below). Additional discussion about the boundary 
conditions for heat transfer problems can be found in Luikov 
(1980) . 

c) Boundary Conditions for Mass Transfer Problems, 
i) Concentration specified at the interface (Dirichlet condition), 

C - C (24). 

a b 

In many instances, the given concentration at the interface is the 
equilibrium value for the system under consideration. 
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ii) Mass flux specified at the surface (Neumann condition). 


( D V C ) ’ n = g(x.,t) (25) 

a 

If there is generation/consumption of solute at the interface due 
to a reaction, g is given by chemical kinetics, 

iii) Mass flux at the surface specified by a mass transfer 
coefficient (Cauchy condition). 


( D V c ) ' n = k ( C - C — ) (26). 

a m a 

iv) Continuity (or discontinuity) of the mass flux at the sur- 
face or phase interface. When solute is generated at a phase 
interface due to a heterogeneous chemical reaction, the mass 
flux has a discontinuity at that interface. One example of this 
kind is the Stefan problem for alloy solidification (Sec(3A.2)). 

In symbols, this condition is then, 


(0>7c) - (D yc) )’(grad F) = (C(F ) - C(F + ))( t) 

a b 


(27) 


Very often, however, Eqn(23) is replaced by Eqn(25), requiring 
only the evaluation of the function g 
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A Note on Flow in Thin Liquid Filins and on Flows with Negligible 
Inertia 

It can be expected that, when a fluid is in the form of a thin 
film, far less time will in general be required for the attainment 
of approximate equilibrium in the direction of the film thickness 
than in the direction of its length. Furthermore, the flow of 
liquids in thin films is an example of a system in which both 
viscous forces and surface tension effects play important roles. 

The governing equations for fluid flow can be considerably 
simplified for the case of flow of thin liquid films. This 
simplification comes about because of the following reasons: 

(i) Since the thickness of the film is small, all velocity 
derivatives across the film are large compared to those along 
its length, and (ii) the flow in thin liquid films can be safely 
assumed to be quasi-unidimensional . 

If a coordinate system is chosen with the x- axis in the 
direction of the length of the film and the y- axis in the 
direction of its thickness, the Navier-Stokes equations become, 

? v /-J t + v 3v /T>X + V 3v /3 y = - (1/P )3 p/^x + 
x xx y x 

+ (/*//> )3 2 v x /'3y 2 + f x (28a) 

and 

'b p/*dy = 0 (28b) 
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If one of the surfaces of the film is a free surface, the pres- 
sure there can be assumed equal to the capillary pressure and if 
the film is so thin that virtually the same pressure exists across 
its thickness, the first term on the R.H.S. of Eqn(28a) becomes. 


- (1/^0 ) dP r /dx = (O'/ ) d 3 H/dx 3 (29) 

Now, since v = ^ H/3t + v 3 h/“Jx , the use of the 

y x 

equation of continuity, for steady state conditions, leads to, 

?H/?x = -(l/v)'B( fv dy)/3x (30) 

X • X 

Equations (28) -(30) constitute the mathematical representation 
of the fluid motion which takes place inside thin liquid films. 

It is well known that the presence of the non linear term 
v'Vv in the Navier-Stokes equations for the thin film makes 
the solution very difficult for all except the simplest flows. 

It is possible, however, to neglect this term in some cases. 
Although the term may not be really zero in these circumstances, 
it will be relatively small and the approximation can be justi- 
fied. Furthermore, if steady state exists, the entire inertia 
term in the equations can be neglected when compared with either 
the pressure or the viscous forces. Flows such as these are 
known as flows with negligible inertia. The geometrical arrange- 
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ment characteristic of flows with negligible inertia is such that 
the thickness of the film must vary along its length and that the 
tendency of the motion must continually be to drag a supply of fluid 
from the thicker to the thinner portions of the film. 

The formulation resulting from the use of the assumptions men- 
tioned above is known as the theory of lubrication. In summary, the 
basic assumptions of the theory (e.g. Hamrock and Dowson(1981) ) are, 

i) the inertia and body forces are negligible compared with the 
pressure and viscous terms, 

ii) there is a negligible variation of the pressure across the film 
thickness, 

iii) the derivatives of the velocity in the direction of the thick- 
ness are much larger than the derivatives with respect to the film 
length, 

iv) laminar flow conditions exist, 

v) the fluid properties are constant across the thickness, and 

vi) there is no slip between fluid and solid boundaries. 

In the most general case, the introduction of the assumptions 
of lubrication theory into Eqn(13) leads to 


v p 



(31) 


The introduction of suitable boundary conditions and the consid- 
eration of Eqn(12) lead to solutions of the general form, 
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(32) 


v/V = v/V (x/L, geometry of the B. C.'s ) 

— o — o 

if the boundary conditions involve v only. Here V q is a 
reference velocity, usually that of the moving substrate. 

For additional information on lubrication theory one can 
consult Batchelor or Schlichting (1979) . 

For the verification of approximations such as those of 
lubrication theory it is frequently necessary to determine the 
relative importance of the various terms in the governing 
equations. For this purpose one can use a so called order of 
magnitude analysis (see e.g. Schlichting). The final result 
of such an analysis is in the form of certain quantities called 
dimensionless numbers. These numbers are actually ratios de- 
scribing the relative importance of the various terms in the 
governing equations. The actual numerical values of these 
ratios can be an useful first guide towards the understanding 
of complex physical processes involving heat transfer and 
fluid flow. Since RSP systems are characterized by these 
features, in Table (1) below we list some of the dimensionless 
numbers more frequently found to apply. For additional 
details on the subject of dimensional analysis the reader can 
consult the book by Szekely (1979) . 
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Table (3A.1.1),- Dimensionless Numbers Useful in RSP 
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heat transfer by conduction 


Table (3A.1.1).- (contd . 
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3A.2.- 


The Formulation and Solution of Solidification Problems. 


Introduction 

The basic feature of solidification problems (also called 
Stefan problems) is that they are represented by a parabolic 
diffusion equation which has to be solved inside a region whose 
boundaries are to be determined as part of the solution. In the 
typical solidification (melting) problem, a substance has a 
transformation temperature at which it changes phase with emission 
or absorption of heat. The liberation or absorption of heat 
takes place at the moving surface of separation between the two 
phases. This surface of separation together with the temperature 
field inside the two phases constitute the solution to the 
solidification (melting) problem. 

Stefan problems are at the core of casting metallurgy and good 
reviews of the metallurgical aspects are available (e.g. Flemings ( 
1974)), Here we will concentrate on the mathematical aspects of the 
subject by presenting the formulation of the governing equations 
and the available techniques for their solution. 

Formulation 

The mathematical statement of the Stefan problem, in the absence 
of fluid motions, consists of the following sets of equations: 
i) The heat equation 
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( 1 ) 


= div( oC grad u ± ) + r h 

i = s , 1 

ii) the heat balance at the moving boundary (Stefan condition) , 

1 

J> L* F/'bt = ( K grad u ) * grad F (2) 

8 

where F (x^ t) = 0 is the phase change surface. 

iii) Equilibrium phase change temperature at the interface 

u = u = u (3), 

s 1 f 

iv) appropriate boundary conditions on all other boundaries, 
which could be of Dirichlet, Neumann or Cauchy type. 

v) Initial conditions. 

Detailed reviews of the formulation of Stefan problems can be 
found in Tayler(1975) and in Crank(1984) . 

The problem represented by (i)-(v) above is nonlinear because 
of the moving boundary. The temperature field depends on the 
exact location of the boundary and this in turn depends on the 
temperature. 

Exact, closed form solutions to the Stefan problem are known 
for only a small number of simple situations. In all cases these 
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solutions apply to one dimensional situations. The one dimensional 
Stefan problem is described by the following equations (see e.g. 
Carslaw and Jaeger (1959) ) , 

'Ju /-$ t = 'b( OC ( “d u /*dx))/'8 x (4a) 

s s s 

? u^t = ^ *3 u^/^x))/^ X (4b) 


K 9u /3 x - K u /3 x = p L 3 X/3 t (4c) 

s s 11 


and 


u 

s 


(4d) 


Closed Form Solutions 

Only two exact solutions of Eqn(4) are available. These are the 
solution due to Neumann and the one due to Schwarz. Both are 
applicable only to infinite regions and under Dirichlet type 
boundary conditions. Schwarz T s solution , however, incorporates the 
mold into the calculation. 

a) Neumann T s solution. If the problem represented by Eqns(4a)-(4d) 
is complemented by boundary and initial conditions such as 

u(0,t) = u 1 < u f (5a) 

and 

u(x,0) = u 2 > u f (5b) 
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the solution is 


u g = + (u^ - u^) erf( x/2(oC g t)* 72 )/ erf A (6a) 

u 1 * u 2 + (u f - u 2 ) erf c ( x/2(e * ^ t) 172 )/ erfc(A( ot/ ^) 1 ^ 2 ) 

(6b) 

X = 2 A ( <* s t) 1/2 (6c) 

where the quantity A is the root of 

^ oC l J 2 (u 2 - u f ) exp( - X 2 / 0^) 

K s OC J /2 (u f - U x ) erfc ( A( oC^^) 172 ) 

- X 1 1rl/2/ C p s <u £ - V (6d) 

In Sec (5 ,2) above we have included a FORTRAN program called 
NEUMANN which computes the instantaneous moving boundary location 
X (t ) and the temperature fields, u (x,t) , and u (x,t) by means 

S 

of Eqns(6a)-(6d) . 

b) Schwarz's solution- Here the chill is incorporated into the 
problem by imposing Dirichlet conditions at a point in the mold 
far away from the mold-casting interface. Figure (1) shows 
schematically how the two solutions differ. 


exp( - A 2 ) /erf A - 
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Since now the mold is included in the calculation, the heat 
equation for the mold must be solved together with Eqns(4). For the 
mold we thus have 


^u m /3t = *3 ( ©C m ( ^u m /^x))/3x (7a) 

with the boundary condition 

u m -*■ U 1 < u f , when x -* - mo (7b). 

The solution of the system of Eqns(4a)-(4d) , (7a) and (7b) is 


u 


K s< /2 <"f - u l> 


m 


K OC 1/2 
s m 


+ K oC 1/2 
m s 


erf /V. 


( 1 + erf ( x/2 ( oC t) 1/2 )) + 
m 


+ u. 


for x < 0 


(8a) 


(u, - 


u 


+ '/f erf( x /2(oc t) 1/2 )) 


K OC 1/2 
s m 


K oC 1/2 erf \ 
m m 


+ u x > for 0 < x < X (8b) 

(u f - u 2 ) erfc(x/2(ot i t) 1 ^ 2 ) 

U 1 = " + u^ , for x > X 

erf c ( A ( OC / 0O 1/2 ) 

S X 

(8c) 

X = 2 OC g t) 1/2 (8d) 
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where A is the root of 


K m*s /2 + 

K 0C 1/2 + K OC 1/2 erf A 

s m m s 

K 1 ^J 2 (u 2 " U f ) ® XP( " *s ^ 2/ ° t l ) L ^ 1/2 A 

K (u f - u ) erfc ( ^. ( C (u -u ) 

slrl si ptl 

s 

(Be) 

In Sec (5. 3) we have included a program called SCHWARZ which 
computes the instantaneous position of the moving boundary and 
the temperature field according to Eqns(8). 

c) The Case of Alloys. While pure substances solidify (melt) at 
a fixed temperature in normal circumstances, alloys change phase 
along a temperature range. This range is limited by the liquidus 
and solidus temperatures of the alloy, respectively, u and 
Ug . Moreover, during alloy solidification more or less severe 
solute redistribution processes take place. Here we restrict our- 
selves to the thermal problem which can be handled in a way 
analogous to Neumann’s and Schwarz's solutions above if the 
following simplifying assumptions are introduced; 

i) the heat of freezing is liberated uniformly over the melting 
range, 

ii) the liquid is initially at the liquidus temperature, and 


180 


iii) inside the melting range a specific heat given by 

C eff ' C Pl + L/ < u l " “s' <«>• 

is used. 

It is now possible to write the solution to the alloy solid- 
ification problem according to either the Neumann or the Schwarz 
formulations. However, instead of Eqns(6d) and (8e) one should 
use 


exp( ( oC g - O^) Xf / OCp erfc ( X ( OCJ oC^) 1 ^ 2 ) 


erf 


< U L - U S } K 1 * s /2 

(u s - u i> “- 1' 2 


( 10 ), 


or 


1/2 


K m exp( ( oC g - ocj X*/ o^) erfc ( X(o£j O^) 172 ) 


K 1/2 + K O C 1/2 erf A 

s m ms 


<u l - u s > 

(u s - u i } K s 01 1 /2 


an 


These expressions were obtained by simply making L = 0 in 
Eqns(6d) and (8e) (since the latent heat is included in C 

ef f ' • 
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For actual computations one uses instead of and u^ 

instead of u^ in Eqns (6a) -(6c) , (10) or (8a) -(8d) , (11) • The 
programs NEUMANN and SCHWARZ given in chapter 5 can be 
used to perform these calculations after the incorporation of 
these changes. 

Additional analytic or semi analytic techniques have been 
used for the solution of Stefan problems in other geometries and 
with different boundary conditions. However, methods such as 
integral profile, series expansions and invariant embedding are 
limited to still relatively simple configurations. The availa- 
bility of digital computers has led to the displacement of all 
these techniques by more convenient finite difference or finite 
element methods. The analytical solutions remain important, 
however, since they are used to verify the accuracy of numerical 
methods and sometimes to start the computational algorithms. 

The solution of most of the solidification problems found in 
practice will invariably require the use of numerical methods. 

Numerical Methods for Solidification Problems 

Since the advent of digital computers many complex solidification 
problems have been solved. Many numerical techniques have been 
proposed to handle the nonlinearities introduced by the moving 
boundary. No single "best" method seems to exist, however. Gener- 
ally speaking all numerical methods subdivide the region of 
interest into small volume elements. Discrete forms of the 
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energy equation are then written for each element and the entire 
set of resulting equations (for the whole system) is solved by 
standard algebraic methods. 

Broadly speaking, the numerical methods used for the solution 
of Stefan problems can be classified according to the following 
three criteria: 

a) According to the treatment of the moving boundary. In this 
case one has; 

i) front tracking methods (Hsu et al(1981)). By careful readjust- 
ment of the computational grid during the computation, the 
precise location of the boundary is recorded for every time step. 
The computational grid itself is redefined at each step on the 
basis of the motion of the solidification interface, 

ii) Fixed domain methods (Voller and Cross (1981) ) , The Stefan condi 
tion (Eqn(2)) is absorbed into the heat equation by introducing the 
enthalpy. The resulting equation is then solved in a fixed grid. 

The boundary can also be fixed, however by performing a coordinate 
transformation using the solidification interface location as the 
reference length of the transformation. 

b) According to the discretization technique employed. Here we have 

i) finite difference methods (White(1983) ) . These are the easiest 
to implement, however, special equation are required at the outer 
boundaries in the case of irregular domains. 

ii) Finite element methods (Ettouney and Brown(1983) ) . Here the 
discretization equations is derived from a variational principle. 
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Since quadrilateral grids are not mandatory when using F.E., it 
is easier to fit the grid to the boundaries of irregular compu- 
tational domains. 

c) According to the procedure used to advance the solution in 
time. In this case we have; 

i) explicit methods (Voller and Shadabi(1984) ) . The temperature 
at a given location depends only on the temperatures of neigh- 
boring points at the previous time step. 

ii) Implicit methods (Elliott and Ockendon(1982) ) . The tempera- 
ture depends on the temperatures of neighbouring points at the 
present time step. 

iii) Semi-implicit methods. A combination of (i) and (ii) above. 

Much more additional information about all these techniques and 
others can be found in the recent monograph by Crank . 

When selecting a numerical method for a given problem, accuracy, 
ease in programming, stability and consistency are among the most 
important considerations. Since for the work reported in this 
thesis we used an explicit finite difference fixed domain method 
for the solidification calculations, this we will review in some 
detail. For the presentation we follow Elliott and Ockendon 
where the reader will find many details not covered in this brief 
review. 

Fixed Domain Methods 


a) Formulation. When solving problems with moving boundaries one 


is usually interested on the properties of the phase chage boundary. 
Since heat is released or absorbed at the phase change temperature, 
this gives rise to discontinuities in the derivative of u across 
the interface. To deal with this lack of differentiability one 
introduces the concept of weak solutions. To describe weak 
solutions we start by noticing that Eqns(l) and (2) above repre- 
senting the differential energy balance can be written in a 
different but equivalent form by using the enthalpy. From 
thermodynamics , the enthalpy is given as a monotonically in- 
creasing function of the temperature, i.e. 


E = f( u ) 


( 12 ) 


It can be shown that the original Stefan problem described by 
Eqns(l) and (2) can be reformulated in terms of Eqns(12) and (13), 

"BE/3t = div(K grad u ) (13). 

The classic solution to the Stefan problem thus being a pair 
of functions of time and position ^u, such that Eqns(12) 

and (13) , hold. The introduction of the enthalpy is useful for 
computational purposes since the location of the phase change 
surface is now implicit in the governing equations. Thus, instead 
of focusing on the pair |u, F j , we concentrate on ^u, Ej , 
with a considerable programming simplification. One can readily 
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verify that the formulation of the Stefan problem in terms of u 
and E is entirely equivalent to the original formulation in 
terms of u and F through the introduction of the concept 
of weak solutions. 

A weak solution of the Stefan problem is defined as a pair of 
bounded integrable functions ^u, inside the domain of 

interest , such that Eqn(12) is satisfied and the integral identi- 
ty (with Dirichlet data) , 


1 1 ^ & /7>t + ~ 2 


u<7 $ ) dx^ 


dt 


[f e o 0l) d 2L + 0 /'d n) dA 


n. 


dt (14) 


where t is the time and the space, holds for all test 
functions with continuous derivatives •9^/a t , aZ/Kx, 

and ^ ^ , such that 0 = 0 on x. = 0, 1 

The two most important properties of such weak solutions are 
the following. Firstly, it can be proved that the weak solution 
exists and is unique. Second, certain difference schemes converge 
to the weak solution. These two properties are an essential 
requirement for any numerical technique to be useful. Most im- 
portantly, the consideration of weak solutions eliminates the first 
spatial derivative of u from the formulation making then 
unnecessary the separate consideration of solid and liquid regions 
in the computer code. The calculations are instead performed inside 
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the entire domain of interest using the same finite difference 
form and the location of the solidification interface is determined 
from the resulting temperature field. So, no complicated front 
tracking procedures are necessary and a rectangular grid can 
conveniently be used. More information about the mathematical 
aspects of weak solutions can be found in Atthey (1975) . 

b) Discretization of the weak formulation. Equation (13) must be 
put in a form suitable for computer calculations. This process is 
called discretization. Discretization can be done according to 
any of several more or less standard procedures. Finite differences 
and finite elements being two popular examples of discretization. 
Since for the calculations reported in this thesis a finite 
difference method was used, this we will describe. The main ideas, 
however, are independent of the discretization method used. 

Assume first that the domain of interest is covered by a 
uniform , rectangular net. A finite difference approximation to 
Eqn(13) can be obtained by simply replacing spatial derivatives 
with central differences and the time derivative with a one 
sided difference (see e.g. Ames (1977) ) . The discrete problem 
thus consists of finding vectors _E n + 1 and u 11 + 1 , the 

values of the enthalpy and the temperature at the mesh points 
for the n + 1 time step such that. 


E 


n + 1 


+ 


A 


n + 0 
u 


0 (15) 
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where 


n + 
u_ 

and such that 


= e u n + 1 + (i - e ) u n 


E 


n + 1 
i 


e , n + 1 . 

f ( u ) 


(16) 


(17). 


Note that, in Eqn(15) the matrix A is the finite difference 
approximation to the Laplacian operator. Also note that the values 
of 0 = 0, 1/2, and 1 , correspond, respectively to the 

explicit, Crank-Nicolson and implicit time discretizations. In 
particular, the explicit scheme has been found to be stable and 
convergent to the weak solution as long as the stability restric- 
tion 

<x A t/( Ax ) 2 £ 1/2 (18) 

is satisfied. The stability constraint, however, may lead to 
prohibitively small time steps when using fine grids. In this 
circumstances the added complexity of an implicit scheme may be 
warranted. Elliott and Ockendon have proposed a successive 
overrelaxation algorithm for the solution of the equations 
resulting from the implicit scheme and they have proved that 
their method is stable and convergent. The main reason for the 
adoption of the explicit scheme in this work was, however, the 
simplicity of its implementation. 
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c) Weak solutions for alloy solidification problems. Most materials 
of metallurgical interest are alloys and it is indeed fortunate that 
the enthalpy method can be readily extended to deal with alloy 
solidification problems. Since alloys do not solidify at a fixed 
temperature but along a melting range, the enthalpy— temperature 
curves for alloys do not show the jump characteristic of pure 
substances (see Fig(2)). The explicit finite difference scheme 
of the weak formulation for the alloy solidification problem is, 


with 


,n + 1 


p n a n 

E — A u 


(19) 


n + 1 


“ <E i + 1 - E f> 


+ u. 


= V 


< 1 ' e ( >(«l-u s > + “< 




n + 1 


u„ 


for E" +1 > E f 
. , 0 i E" +1 i E f 

Ef 1 < 0 


( 20 ) 

In this set of equations, however, one usually uses u = K u 

i i 

The computational procedure is as follows: Eqn(19) is used 

first to calculate E n+1 in the entire domain, then, Eqn(20) 
permits the calculation of the temperature field u n+ ^ . This 
algorithm is stable and converges to the weak solution of the 
original problem as long as the stability constraint given by 
Eqn(18) is satisfied. 

It must be said, however, that the alloy solidification problem 
also requires the calculation of the solute distribution resulting 
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from the difference in solubility of the impurity in the solid and 
liquid phases. Thus, a species transport problem has to be solved 
simultaneously with the heat equation. So besides solving Eqns(19) 
and (20) one also has to solve (see Sec(3A.l) and also Wilson et al 
(1984)), 

^ C./? t + (^ C )' v = div( D grad C.) + r (21) 

1 i l c 

i = s , 1 

where is the concentration of solute in the solid or 

liquid phases for the particular case of a binary alloy system. 
Moreover, a "jump" condition, analogous to the Stefan condition 
for heat transfer, must also be considered in the solution of the 
mass transfer problem, i.e., 

( C(F ) - C (F ) K^F/^t) = ( D grad C ) * grad F 

1 

( 22 ). 

Furthermore, since the heat and mass transfer problems are 
coupled, one usually assumes thermodynamic equilibrium at the 
solidification interface, i.e., 

u(t,F") = u(t,F + ) = m g C(t,F~) + u A 

= C(t,F + ) + u A (23) 

where u A is the melting point of the pure solvent. 
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Wilson et al. have reformulated the alloy solidification problem 
in terms of weak solutions in an analogous way to the enthalpy 
formulation of the Stefan problem. They claim good results from 
the use of an explicit finite difference scheme to solve the 
coupled heat and mass transfer problem. 

d) Solidification in the presence of fluid motion. This problem 
has been the subject of continued research for at least the last 
30 years. The peculiar effects of fluid flow on the solidified 
material were recognized early on (e.g. Flemings(1956) and Roth and 
Schippers(1956) ) . The role of fluid motion in the dissipation of 
the melt superheat was also noticed (Adenis et al(1962)). Sahm(1982) 
has presented a comprehensive survey of the metallurgical effects 
of fluid motion during solidification. 

The mathematical problem of phase change with fluid motion 
involves the solution of the fluid flow equations together with 
the energy balance incorporating a convective term. This is not a 
straightforward problem since the indetermination in the location 
of the solidification interface makes the flow problem into a 
non linear boundary value problem with unspecified boundaries. 

Thus, although the existence of a weak solution to the Stefan 
problem with convection due to Stokes' flow has been mathemati- 
cally proved (Cannon et al(1983)), the actual numerical solution 
of specific problems is still the subject of active research. 

Most treatments to date have circounvented the flow problem. 
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Typically, a heat transfer coefficient is introduced at the phase 
change surface to account for the fluid motion (Larreq et al(1978)). 
Among the few reports of rigurous solutions most apply to 
simple geometries. Sparrow et al(1977) used a front tracking method 
to model melting inside a narrow cavity with natural convection. 
White (1982) used a fixed domain method to solve a somewhat similar 
problem, 0 T Neill (1983) and Argyris et al(1984) have used finite 
element analysis to model melting with fluid flow in cavities. 

Oreper and Szekely (1984) used a hopscotch fixed domain algorithm 
to describe melting and electromagnetic stirring inside welding 
pools. 

Since the additional non linearities introduced by the flow 
equations will most likely complicate the solution of the Stefan 
problem, simplified versions of the general problem can be 
expected to be more easily solved. For example, the calculations 
by Miyazawa and Szekely (1979) , (1981) and the ones reported in 
this thesis for the Stefan problem with fluid flow produced by 
an inertialess fluid, are a proof of this claim. 



3A.3.- The Solution to the Fluid Flow Equations for the PFMS 
System. 

For the schematic PFMS puddle shown in Fig (3. 3. 1.2) ( and 
photographed in Plate (3 . 3 .1 .1) ) , the governing equations for 
fluid flow at steady state are; 
the equation of motion 



One can integrate Eqn(l) with respect to y for a fixed x— 
location to obtain an expression for , i.e. 


V 

x 


Z J*- ) (dP/dx) y 2 + K x y + K 2 


Where the coefficients and depend on the particular 

boundary conditions used. 

The boundary conditions on the top surface of the puddle change 
after the puddle detaches itself from the nozzle lip in the 
downstream direction. Because of this, the regions before and 
after the detachment point must be treated separately. 


* 


<• 
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e) Points before detachment. In this region the appropriate 
boundary conditions are: 

i) the no— slip condition at the puddle nozzle interface, i.e. , 

V x " v y = 0 at y = H (4a) 

and (ii) the continuity of velocity at the solidification 
interface, i.e., 


V 

x 


v V = V on y = y (4b) 

» y I„ s 


where Vj is the (Eulerian) velocity of the interface. 

Substitution of Eqns(4) into (3) and rearrangement leads to 


where 


A 

o 


A 


1 


and 


A 


2 


v x = A 0 + A i y + a 2 y 2 (5a) 

(1/2 yx ) (dP/dx) H y g - Vj /((y /H) - 1) (5b) 

/ x 

/ (y g “ H) - (1/2 /-*• ) (dP/dx) (y + H) (5c) 

x / s 

(1/2 ) (dP/dx) (5d) 
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We can readily see that the specification of five quantities, 

v , H , ^ , dP/dx , and allows the calculation of 

the velocity field. V_ ( = V ) and H (- H ) are 

I r o 

x x 

usually input data whereas y^ can be calculated from any 
of several solidification heat transfer models. dP/dx , on the 
other hand, has to be computed such that it satisfies the 
equation of continuity . At any given downstream location x, 
the total mass flow rate crossing through the entire puddle 
thickness is given by 


Q = 



( 6 ) 


where Q and Q are, respectively , the flow rates 

carried by the partially solidified ribbon and by the fluid film 

above it. Now, since the ribbon moves like a rigid body while 

the film is a Newtonian fluid, Q and CL are given by, 

s X 


and 



(7a) 


(7b) 


Performing the integrations of Eqn(7) and combining the 
result with Eqns(5) and (6) gives. 
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Q/w = V r y g + V r /2(H - y ) + 

x x s 

- (1/2 AA.) (dP/ dx) (H 3 / 6 ) ( 1 - 3 (y /H) + 3(y /H) 2 - (y /H) 3 

/ s s s 


whence 


dP/dx 


(V /2(H - y )) + V y - q/w 

1 o L S 

12 yu. — 5 

H 3 (1 - 3(y g /H) + 3(y g /H) 2 - (y s /H) 3 ) 


( 8 ) 


(9) 


As expected, in the absence of solidification, Eqn(9) reduces 
to the well known one dimensional form of Reynolds' equation of 
lubrication theory (see e.g. Szekely(1979) ,p.H5) . 

Equations (5) and (9) allow us to compute the velocity field 
once the process parameters and the solidified thickness are 
known. 

Sometimes it is useful to present the results of flow calcu- 
lations in terms of the stream function rather than the velocity. 
A normalized stream function can be defined for our system as, 



( 10 ) 


Substitution of Eqn(5) into (10) gives. 


r 


= w ( V 


r Y s + (A 2 /3)(y 3 - y 3 ) + (A /2)(y 2 - y 2 ) + 

X X S 


+ A o (y - y s ) )/ Q 


(ID 
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b) Points after detachment. After detachment the no slip condi- 
tion can not be applied to the top surface of the puddle. This 
is so because the fluid particles on this melt-gas interface are 
much more mobile than those on the melt-nozzle interface. For 
this case the proper boundary conditions are then, (see Sec(3A.l)) 

dV^/dy = 0 on y = H (12a) 

and 

V = V = V 

x I x r x on y = y g (12b) 

Equation (12a) simply represents the continuity of shear force 
across the melt-gas interface while Eqn(12b) is again the no 
slip condition applied to the solidification interface. 

It must be noted that the precise location of the melt-gas 
interface (the free surface) is not known in advance but must 
be calculated simultaneously with the velocity. The calculation 
of H after detachment requires the consideration of the 
capillary effects due to the surface tension of the melt and the 
details are described in Sec(c) below. For the time being, we 
shall assume that such calculation has already been performed 
and that the function H(x) is known. Under these circumstances, 
the combination of Eqns(12) with Eqn(3) gives, 

v x = A q + A x y + A 2 y 2 (13a) 
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where 


A = 
0 

V - (1/2/x ) (dP/dx) (y 2 - 2 y H) 

r x / S s 

(13b) 

A 1 = 

- (l/2ykA) (dP/dx) (2H) 

(13c) 

A 2 " 

(1/2 M ) (dP/dx) 

(13d) 


Apart from the fact that the value of H(x) must be obtained 
from a separate calculation, Eqns(13) are analogous to Eqn(5) 
for points upstream. 

Since Eqns(6) and (7) are still valid after detachment, they 
can now be combined with Eqn(13) to obtain. 


Q/w = V y + V (H - y ) - (1/2/*.) (dP/dx) (2H 3 /3) * 

t 5 L S / 

XX ' 

*( 1 - 3(y /H) + 3(y /H) 2 - (y /h) 3 ) (14) 

o s s 

and 

V r (H - y g ) + V r y s - Q/w 

dP/dx = 3M — ^ ^ ( 15 ) 

H 3 ( 1 - 3(y /H) + 3(y/H) 2 - (y /H) 3 ) 

o S S 

The corresponding expression for the stream function can now 
be derived. This turns out to be identical to Eqn(ll) except that 
the coefficients A^, A^, and A£ are now given by Eqns(13). 
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c) Calculation of the meniscus shape. The region around the free 
surface can be considered as constituted by three parts, namely 
(i) the melt, (ii) the gas phase, and (iii) the interface between 
them. In such a system the changes in free energy are related to 
the changes in the volume of the bulk phases as well as to the 
changes in the surface area of the interface. Denoting the free 
energy change by dF , we have that 


dF 



- P ) dA dH 1 



(d 2 H 1 /dx 2 ) dA dH 


(16) 


Here, P^ and are the pressures inside the bulk phases and 

dA is an element of area of the interface. For equilibrium, the 
free energy is a minimum, i.e, dF = 0 , thus 


P 1 " P 2 = ' O' ( d2 V dx2 > 


(17) 


The difference in pressure between two contiguous phases 
separated by an interface with surface tension <r is an important 
quantity in capillary hydrodynamics called the capillary pressure 


= - O' (d 2 H 1 /dx 2 ) 


(18) 


As described in Sec(3A.l), the equation of motion for the thin 
film of melt contained between the meniscus and the solidification 
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interface, under the assutnpt ions of lubrication theory, is 


dP/dx ~ M d 2 V /dy 2 (19) 

/ X 

If the liquid layer is so thin that the pressure is essentially 
constant across its thickness, the substitution of Eqn(18) into 
Eqn(19) leads to. 


(S~ (d 3 H 1 /dx 3 ) 


+ JJ* (d 2 V x /dy 2 ) = 0 


( 20 ) 


Equation(20) can now be integrated with respect to y to 
obtain a cl®**d form expression for . This expression, in 

turn, can be combined with Eqns(6) and (7) above to produce an 
equation for the fluid film thickness . This is 


d 3 M 1 /*c 3 - (3 M /<T )((Q 1 /wH 3 ) - V r /H 2 ) = 0 (21) 

Recall that here, as before, . = H - y^ * f° n (2 1 ) is 

a third order non linear ordinary differential equation which can 
be used to compute the precise location of the meniscus. Closed 
form solution to Eqn(21) are known only for the asymptotic cases. 
Therefore, a w*iserical method is required to solve it for the 
conditions of our system. This we describe next. 
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To solve Eqn(21) using numerical methods we first transform it 
into an equivalent set of three first order equations, i.e., 

d^/dx = Hj 

dH j/dx = H” 

and 

dH"/dx = (3M/(T)((Q 1 /whJ) - V r /H* ) (22c) 

' X 

This system can be written in the abbreviated form 

dH^/dx = f_(x) (22d) 

where H_^ = ( H^, Hj, ) and f_ is the vector formed by 

the right hand sides of Eqns(22) . 

Although high order schemes may prove to be more accurate, for 

the calculations reported in this thesis we have used the simplest 

approximation method for initial value problems, namely Euler's 

forward method. As initial conditions for Eqn(22) we have the 

values of , Hj , and at the detachment point. The value 

of H_ = H - y , depends on the solidified thickness at the 
1 o J s 

point of detachment. The values of and H” , on the other 

hand, have been estimated from still frames of high speed movies 
of the puddle during PFMS . Our estimates, from the photos, 
were , 


(22a) 


(22b) 
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= - tan Q (23a) 

and 

" 2/ L p (23b) 

The discrete analogue of Eqn(22d) , obtained using Euler's method, 
is (see e.g. Dahlqvist and Bjtfrck(1974)) 

-1 +1 = -1 + Ax l (x) (24) 

The calculation is then advanced step by step in the downstream 
direction. Starting from the values given by Eqns(23) , the 
repeated application of Eqn(24) allows the calculation of the 
location of the free surface. However, note that the solidified 
thickness must be calculated prior to solving Eqn(24) at each 
step to account for the solidification. It is here where the 
coupling between the flow problem and the solidification problem 
is made (see Sec (3. 3.1)). 
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4A,1.- A Comment on a Numerical Method for the Solution of Problems 
Involving Transport Phenomena. 

Introduction 

When facing a problem involving transport processes it is always 
advisable to look for suitable simplifying assumptions capable of 
reducing the mathematical complexity of the problem. In many cases, 
however, this is not possible and one must resort to numerical 
methods for the solution of the governing partial differential 
equations. The advent of powerful computers has contributed a great 
deal to the development of the new field of numerical heat, mass, 
and momentum transfer. Many procedures have been proposed to deal 
with the equations of transport phenomena. Unfortunately, however, 
only a few of them are in the form of commercially available, general 
purpose computer programs. One of the most successful methods is 
the one developed by Patankar and Spalding, among others, at the 
Imperial College of London. The Patankar-Spalding (P-S) method solves 
the equations of transport using implicit finite difference schemes 
derived from a control volume formulation. In the next few pages we 
present a brief description of some of the main features of this 
method hoping to familiarize prospective users of the commercially 
available version. For additional information the reader can 
consult the presentations by Patankar (1980) and Spalding (1980) . 


The General Transport Equation 


o 


A careful look to the transport equations presented in Sec(3A.l) 
reveal that they all have very much the same mathematical form 
So, if instead of the physical variables velocity, temperature, and 
concentration we introduce the generalized transported variable 


* 


, the differential balance, expressing the conservation of 
can be written as 


v = div( r ) + 

where I"* is the diffusion coefficient for * and 
the source term accounting for any absorption-release of 


S (1) 



inside the system. 

The recognition of the mathematical similarity of the various 
conservation equations produces considerable simplification of the 
computational procedure since more or less the same method can be 
used to find v_ , u, and C . It must be also recognized, 
however, that the solution of the fluid flow problem is more in- 
volved than that of the heat and mass transfer problems. The added 
complexity comes about because of three main reasons. First, the 
Navier-Stokes equations (Sec(3A.l)) are non linear. Second, four 
equations must be solved to determine the flow field, 3 for the 


velocity components and one for the pressure (continuity). Finally, 
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special techniques are required during the solution of the flow 
equations to be ablo to obtain physically realistic solutions. 

It is indeed fortunate that despite these complexities the P-S 
algorithm retains its simplicity. 

The basic idea of the numerical method is the replacement of 
the governing equations by simple algebraic analogues which can 
in turn be solved using a digital computer. The two main tasks are 
then, first to devise a method for the derivation of the discrete 
equations, and second to find a convenient method for their solution. 

Discretization of the Governing Equations 

To discretize the governing equations one starts by mentally 
subdividing the region of interest into a number of small domains 
called control volumes. The governing equations are then integrated 
over the extent of such control volumes. Finally, one assumes a 
local linear variation of the field variables inside every volume 
element to arrive at the final set of algebraic equations repre- 
senting the conservation principles inside the control volumes in 
a discrete sense. The solution of this set of algebraic equations 
is the (approximate) solution to the original problem. 

Without going into details and for the sake of illustration we 
now present the typical form of the resulting algebraic equations 
for the control volume centered in P shown in Fig(l), as given 
by Patankar for the case of unsteady two dimensional heat con- 
duction with internal heat generation. 
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( 2 ) 


a P U P 


a E U E 


a w u w 


+ a N U N 


a s u s 


4- b 


where 

a E - K e &y/(Jx) e , a y ' K „ Ay/<;«)„ 

a N ‘ ' a S ' K S 4l,( h> s 

a° = P C A x Ay/ £t 
P * p 

b = S c A x Ay + a° u° 

and 

3 P = a E + a W + a N + a S + a P “ S P Ax Ay 


(2a, b) 


(2c, d) 


(2e) 


(2f ) 

(2g) 


In Eqns(2), Ax Ay *1 is the volume of the control volume 
and A t the time step . Moreover, the source term has been 
linearized, i.e. S = u p • Note that if the domain 

has been divided into N control volumes, the problem at this 
stage consists of solving a system of N algebraic equations 
with N unknowns ( i=l,...,N). 

The Effect of Fluid Motion 

When the field variable being transported travels by convection 
as well as by diffusion the discretization equations must be 
slightly modified to account for the fact that convective transport 
takes place basically in the downstream direction of the flow . 

This procedure is called upwinding . 
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To deal with upwinding, Patankar proposes the so called "power 
law" scheme based on the results of comparison of the results 
obtained from it with the closed form solution of the one dimen- 
sional convection-diffusion equation. It is to the credit of the 
P’-S method that the discretization equations remain basically 
unchanged after the introduction of upwinding, except for some 
changes in the actual values of the coefficients. So, for the 
two dimensional convection-diffusion problem the actual coefficients 


are; (see also Fig(2)) 

a E - D e A(|Pe e l) + [ [ - F g , 0 ] ] (3a) 

a w = D w AdPeJ) + [[ F w , 0 ] ] (3b) 

a N = D n A(|Pe n |) + [[ - F n , 0 ]] (3c) 

a s = D g A(|Pe s l) + [[ F g , 0 ] ] (3d) 

a° = C p Ax Ay/ At (3e) 

b = S c A x Ay + a° u° (3f ) 

a p = a E + a w + a N + a s + a p “ s p Ax Og) 


where 
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(3h,i) 


D e = K e Ay/(5x) , D « K Ly/(Sx) 

e e e ww w 

D n = K n Ax/( ^ y) n * D s = K s Ax/( ^ y > s Oj.k) 

A(iPe|) = [[ 0,(1- 0.l|Pe| ) 5 ]] (31) 


Pe e = F e /D e ’ Pe . 


w = F „ /D „ » Pe „ = FVd , Pe = F /D 

w ww n nn s ss 

(3m,n,o,p) 


and 


F = (O v ) A y 
e j x e 


F w ■ \>W Ay 


(j»Vn Ax 1 F s = 


(3s , t) 


In Eqns(3) , the symbol [[ , ]] is a well known FORTRAN 
operation which simply selects the greater of the two quantities 
separated by the comma . 

Note that even now the problem still consists of solving a system 
of as many equations as volume elements there are in the domain. 


The Computation of the Flow Field 

Many of the difficulties associated with solving the flow 
equations have been solved by the use of staggered grids of the 
type shown in Fig(3). The discrete equation for the pressure is 
derived by combining the discrete versions of the equations of 
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continuity and motion. The resulting set of equations for the 
general heat, mass, and momentum transfer problem consists of the 
the following (in the case of N control volumes): (a) One set 

of N equations for the pressure field, (b) three sets of N 
equations each for the three components of velocity, (c) one set 
of N equations for the temperature field, and (d) as many sets 
of N equations each as chemical components there are in the 
system. In practice, many transport processes take place under 
turbulent conditions and in this case additional conservation 
equations (for turbulence quantities) must be introduced (from 
corresponding turbulence models) to account for this. The references 
given should be consulted for details in this regard. 

The Solution of Systems of Algebraic Equations 

From all the above it should be clear that once the discretiza- 
tion procedure is finished the problem has been reduced to solving 
sets of algebraic equations. These equations are not necessarily 
linear. The method of solution, however, is the same regardless 
of this fact. The method is a combination of direct and iterative 
techniques and it is known as the Gauss-Seidel line by line 
method . 

In the Gauss-Seidel line by line method one starts by selecting 
a direction for sweeping. Then , the equations for the grid points 
lying on a line perpendicular to the sweeping direction are solved 
simultaneously using a direct method (e.g. Gauss elimination). The 
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values used for the field variables of grid points on neighbouring 
lines are the latest ones available in the computer memory. Once 
this is done one moves forward to the next line of grid points 
along the sweeping direction and does the same thing. This opera- 
tion is continued until the entire domain has been swept line by 
line. This concludes one sweep. However, since one uses guessed 
values for the field variables in order to start the computation, 
the sweeping operation must be repeated until the values of the 
variables stop changing appreciably from one sweep to the next. 
These final, converged values are considered the numerical solution 
of the original problem. In Sec (5. 7) we have included a program 
constructed based on the P-S method and which is capable of 
computing heat transfer by conduction-convection for two dimen- 
sional situations . The flow field in this case, however, is not 
calculated numerically but from closed form expressions given by 
lubrication theory. Moreover, in Sec (5. 6) we present a very simple 
program capable of solving systems of algebraic equations directly. 
This program can be used as the core of an algorithm based on the 
P-S method. 

In practice, the need to solve many sets of equations simulta- 
neously requires that some thought be put on the sequencing that 
must be followed. Patankar proposes the following series of steps 
as a general algorithm for the numerical solution of problems in 
fluid flow and heat and mass transfer: 
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a) Guess the pressure field, 

b) solve the discrete Navier-Stokes equations to obtain a first 
estimate of the velocity field, 

c) correct the pressure by using the equation of continuity, 

d) compute again the velocities but this time use the corrected 
pressure, 

e) solve the discretization equations for all other field 
variables using the procedure described before. At the end of 
this step, one sweep for all field variables has been completed. 

f) Steps (b)-(e) are repeated again and again until the values 
of all the field variables in all the control volumes stop 
changing significantly from one cycle to the next. 

Conclusion 

The algorithm described above has been used extensively for the 
solution of transport phenomena problems in mechanical, chemical 
civil and metallurgical engineering. It offers the possibility of 
exploring the behavior of complicated systems with a relatively 
modest amount of effort. Instead of constructing his/her own 
program following the ideas presented above and described in 
detail in Patankar, the reader may choose to use the more convenient, 
commercially available versions. An improved understanding of 
RSP systems can be expected from the use of these ideas. 
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